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— — * Th i s  thesis  describes  the  development,  of  a  method  tor  calculat¬ 
ing  room  heating  and  cooling  loads  using  frequency  response  methods 
to  calculate  heat,  flow  through  mul ti - 1  aycred  room  walls  and  roofs.  - 


Existing  methods  tor  calculating  room  loads  are  discussed  and 
the  nred  to  consider  the  effects  of  thermal  mass  arc  emphasized.  A 
detailed  derivation  of  the  solution  ot  the  transient  one-dimensional 
heat  conduction  equation  for  multi-layered  slabs  which  leads  to  the  - 
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Block  20  continued. 

■J  calculation  of  response  factors  and  conduction  transfer  functions  is 
presented.  Next,  the  special  case  of  a  fixed  temperature  on  one  side 
of  a  multi-layered  slab  and  a  sinusoidal  temperature  variation  on  ihe 
other  side  is  examined  and  the  simplifications  which  result  are 
described.  This  frequency  response  method  for  calculating  heat  flow 
can  be  applied  to  buildings  if  the  outdoor  climate  can  be  represented 
by  a  reasonably  small  number  of  sinusoidal  components.  - 

To  determine  the  extent  of  sinusoidal  variation  in  climate,  time 
series  analysis  techniques  were  applied  to  hourly  measured  weather 
data  from  four  climatically  different  sites  in  order  to  separate  the 
deterministic  periodic  behavior  in  the  data  from  its  stochastic 
behavior.  A  combined  deterministic  plus  stochastic  model  was  found 
to  be  adequate  for  the  climate  variables  studied.  This  model 

included  as  its  deterministic  component  sinusoids  with  frequencies  I 

corresponding  to  the  annual  cycle,  the  diurnal  cycle,  and  the  first  i 

harmonic  of  the  diurnal  cycle.  '  j 

The  effect  on  building  energy  calculations  of  using  only  the  | 

sinusoidal  components  of  climate  was  determined  by  comparing  the 
results  of  energy  calculations  which  used  only  sinusoidally  varying 

climate  to  the  results  obtained  when  hourly  measured  climate  data  was  ; 

used.  Very  good  agreement  was  found. 

Since  ignoring  tie  non-si nusoidal  components  in  climate  varia¬ 
tion  produced  little  error  in  calculating  building  energy  use,  fre-  1 

quency  response  methods  were  judgeH  to  be  adequate  for  calculating 
heat  flow  through  walls  and  roofsr"*The  use  of  this  method  as  part  of 
a  simplified  room  load-calculating  procedure  was  compared  to  results 
from  •  more  detailed  approach,  and  again,  very  good  agreement  was 
f ound  . 

The  use  of  this  simplified  method  tor  calculating  room  heating 
and  cooling  loads  will  allow  building  energy  analysis  procedures  to 
be  implemented  on  microcomputers,  which  may  enlarge  the  number  of 
engineers  who  can  use:  such  procedures  to  improve  the  energy  effi¬ 
ciency  of  their  designs. 
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1  INTRODUCTION 


Background 

In  the  past,  designers  of  heating  and  cooling  systems  for  build¬ 
ings  have  focused  on  two  main  objectives:  (1)  insuring  comfort  in 
the  air-conditioned  space,  and  (2)  minimizing  the  first  cost  oi 
building  mechanical  systems.  A  third,  often  unstated,  objective  was 
(and  is)  that  the  designer  earn  a  respectable  wage  for  performing  the 
design  task.  Furthermore,  the  designer  of  building  mechanical  sys¬ 
tems  has  usually  worked  under  the  constraint  that  his/her  work  must 
not  influence  or  affect  the  building's  architecture.  With  these  com¬ 
peting  objectives  and  constraints,  the  detailed  analysis  of  heat  flow 
and  thermodynamic s  in  buildings  was  economically  impractical.  His¬ 
torically,  the  design  of  heating  and  cooling  systems  for  buildings 
has  relied  on  the  most  simple  calculations.  Designers  typically 
applied  simple,  steady-state  heat  transfer  equations  (sometimes 
adjusted  by  the  use  of  pretabulated  correction  factors)  to  calculate 
heat  gains  and  losses  in  buildings.  Heating  and  air-conditioning 
equipment  was  selected  with  sufficient  capacity  to  maintain  comfort 
under  peak  heating  and/or  cooling  conditions.  Suitable  control  sys¬ 
tems  were  designed  to  insure  that  comfort  would  be  maintained  during 
both  peak  and  off-peak  conditions.  Calculation  of  building  heating 
or  cooling  requirements  or  the  expected  energy  consumption  of  the 
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mechanical  system  under  off-peak  or  part-load  conditions  was  usually 
not  part  of  the  design  analysis. 

Recent  increases  in  fuel  and  utility  costs,  along  with  new 
Government  energy  conservation  regulations,  have  caused  a  shift  in 
building  energy  system  design  objectives  and  have  made  building 
planners  more  aware  of  building  thermal  performance.  The  effects  of 
architecture  on  energy  use  are  being  considered,  and  both  the  annual 
performance  and  the  peak  load  performance  of  building  mechanical  sys¬ 
tems  are  being  evaluated  during  building  design.  These  more  detailed 
analyses  require  a  much  larger  number  of  substantially  more  compli¬ 
cated  cal  cul  atioris.  Computers,  once  used  in  the  building  energy  sys¬ 
tem  design  industry  only  as  research  tools  or  to  produce  tables  of 
correction  factors  for  simple  calculation  methods,  are  gradually 
being  used  in  the  analysis  and  design  of  individual  buildings. 

Detailed  energy  analyses  of  buildings  have  revealed  some  common 
trends.  For  example,  building  energy  systems  which  are  "just  big 
enough"  are  often  much  more  efficient  than  systems  which  are  "plenty 
big  enough."  Also,  estimated  annual  fuel  costs  are  strongly  influ¬ 
enced  by  thermal  mass  in  buildings  and  by  part-load  system  perfor¬ 
mance. 

With  the  increased  interest  in  a  more  detailed  thermal  analysis 
of  buildings  and  their  energy  systems,  a  number  of  building  energy 
analysis  computer  programs  have  been  developed  by  Government 
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I aboratories ,  equipment  manufacturers ,  engineering  associations,  and 
private  engineering  firms.  The  following  two  sections  discuss  the 
methods  used  in  existing  energy  analysis  programs  and  a  proposed  new 
method  for  accomplishing  part  of  the  energy  analysis  task. 

Existing  Building  Energy  Analysis  Procedures 

Most  modern  building  energy  analysis  computer  programs  perform 
hourly  calculations  for  a  1-ycar  period,  using  recorded  hourly 
weather  data  for  a  particular  site.  These  calculations  are  divideo 
into  three  major  parts:  (1)  space  thermal  load  prediction,  (2)  air 
distribution  system  simulation,  and  (3)  central  plant  (boiler/ 
chiller)  simulation. 

Space  load  prediction  involves  calculating  the  amount  of  heat 
which  must  be  added  to  or  removed  from  the  space  to  maintain  a 
desired  space  temperature.  The  effects  of  climate,  building  con¬ 
struction,  space  temperature  control  strategy,  occupancy,  and  Mglit- 
ing  use  patterns  are  some  of  the  factors  considered  in  calculating 
space  loads.  The  need  to  evaluate  the  transient  response  of  the 
space  to  these  1 oad-inf 1 uenci ng  factors  makes  space  load  prediction 
compl i cated . 

Once  space  loads  are  calculated,  it  is  necessary  to  simulate  the 
performance  of  the  air  distribution  system  which  will  supply  or 
remove  heat  from  the  space.  This  phase  of  the  building  energy 
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analysis  calculation  process  usually  includes  steady-state  component 
and  control  system  simulation  based  on  basic  principles  of  thermo¬ 
dynamics  and  heat  transfer.  The  results  from  this  part  of  the  calcu¬ 
lation  are  the  hourly  demand  for  hot  and  chilled  water  and  electri¬ 
city. 


The  final  phase  of  the  building  energy  analysis  process  is  simu¬ 
lating  the  central  plant  components  which  supply  energy  to  the  build¬ 
ing  air  distribution  system.  Components  such  as  boilers  and  chillers 
are  usually  simulated  using  polynomial  expressions;  these  expres¬ 
sions,  which  are  derived  from  manufacturers '  data,  empirically  relate 
component  energy  inputs  to  energy  outputs. 

Note  that  for  the  purpose  of  estimating  ex  |  ec  f  i  o  t  r>i  rgy 

consumption,  transient,  thermodynamic  analysis  is  nensmiy  only  *hon 
calculating  space  loads,  lhis  is  btcaust  t  hi  uf’tr  itv-  <.<  tlemin’s 
of  buildings  (e.g.,  concrete  floors  and  brni  walls;  an  <ml>  hi 
steady  state.  For  example,  solar  energy  absorbed  by  a  building  exte¬ 
rior  wall  must  be  conducted  through  the  wall  and  curved  lvely 
transferred  to  the  room  air  before  it  becomes  a  load  on  the  hu.’o  ng 
air-conditioning  system  (or  before  it  offsets  a  space  heating  load). 

On  the  other  hand,  air  distribution  systems  can  be  simulated  on 
an  hourly  basis,  using  steady-state  assumptions,  since  (1)  the  energy 
storage  capacity  of  the  circulated  air  stream  is  very  small  compared 
to  the  energy  being  transferred  to  the  air  by  the  mechanical  system, 


and  (2)  since  the  response  of  the  system  controls  is  much  faster  than 
1  hour.  Similar  reasoning  permits  the  steady-state  simulation  of 
central  plant  components. 

Because  the  transient  analysis  of  heat  flow  in  buildings  adds 
considerable  complexity  and  computat ional  expense  to  the  building 
energy  analysis  process,  developing  improved  procedures  for  perform¬ 
ing  transient  heat  flow  calculations  is  a  worthwhile  goal.  The 
development  of  an  improved  procedure  is  the  main  subject  of  this 
thesis. 


Existing  Methods  for  Space  Load  Prediction 


The  heat  transfer  problem  in  buildings  centers  around  the  solu¬ 
tion  of  the  heat  conduction  equation.  Almost  all  computational 
schemes  assume  one-dimensional  heat  transfer  for  which: 


3‘~T(x,t)  _  J_  oT(x,t) 
^ y2  a  3t 


[Eq  1-1] 


where  T  is  temperature,  x  is  the  space  dimension,  t  is  time,  and  «  is 
the  thermal  diffusivity.  Given  that  this  differential  equation  can 
be  solved,  a  conductive,  radiative,  and  convective  heat  balance  can 
be  computed  for  interior  and  exterior  surfaces  and  for  the  room  air 
volume.  In  addition  to  the  effects  of  climate,  this  heat  balance 
must  include  energy  inputs  from  people,  lights,  and  equipment,  and 
from  the  space  heating  and  cooling  system. 


Existing  methods  for  predicting  space  loads  differ  in  their 
approach  to  solving  the  heat  conduction  equation.  Gupta,  et  al  .  (ref 
1),  have  conveniently  divided  these  methods  into  three  general 
classes:  (1)  numerical  methods,  ( 2 )  harmonic  methods,  and  (3) 
response  factor  methods. 

Numerical  methods  use  lumped  parameter  approximations  to  the 
heat  conduction  equation  and  were  originally  implemented  using 
resistor-capacitor  circuits  on  analog  computers;  however,  digital 
computer  methods  now  predominate,  and  both  space  and  time  derivatives 
are  approximated  as  finite  differences.  As  is  always  the  case  with 
finite  difference  techniques,  accuracy,  cost,  and  model  stability  are 
all  functions  of  number  of  nodes,  the  time  step  used,  and  the  solu¬ 
tion  method  chosen.  Numerical  techniques  have  the  advantages  of 
being  conceptually  simple  and  able  to  hand 1 e  both  linear  or  non¬ 
linear  boundary  conditions. 

Harmonic  methods  can  be  used  to  solve  the  heat  conduction  equa¬ 
tion  if  the  boundary  conditions  are  represented  as  periodic  func¬ 
tions.  (Climate  can  be  reasonably  approximated  by  a  limited  number 
of  coefficients  of  a  Fourier  series;  see  ref  2.)  These  methods 
require  that  the  building  heat  transfer  parameters,  including  convec¬ 
tion  coefficients,  be  constant  with  time  and  that  radiant  heat 


transfer  be  linearized. 


Mackey  and  Wright  (rets  3,4)  tirsl  used  harmonic  methods  to  cal¬ 
culate  the  heat  gains  through  building  walls  and  through  roofs  which 
separated  a  room  having  a  fixed  constant  temperature  from  a 
sinusoidally  varying  outdoor  temperature.  Walls  and  roofs  were 
characterized  by  a  decrement  factor  and  a  lag  factor  for  the  diurnal 
frequency  arid  its  harmonics.  (Tire  'Method  was  intended  tor  design  aay 
calculations  only.)  Mackey  and  Wright  (ref  5)  also  devised  the  so- 
called  sol-air  temperature,  which  allowed  solar  energy  impinging  ori 
opaque  surfaces  to  be  accounted  for  as  an  equivalent  convective  flux. 
Van  Gorcum  (ref  6),  Muncey  and  Spencer  (refs  7,8),  Pipes  (ref  9), 
Gupta  (ref  10),  Sonderegger  (ref  11),  and  others  have  contributed  to 
the  enhancement  of  harmonic  methods  for  space  load  prediction. 

Response  factor  methods,  which  represent  yet  a  third  approach  to 
solving  t.hc  heat  conduction  equation,  are  in  widespread  use  in  t hi. 
United  Slates  and  Canada.  The  major  advantages  of  these  methods  are 
that  they  are  not  numerical  in  the  sense  of  finite  differences  tech¬ 
niques,  and  they  do  not  require  that  the  heat  conduction  boundary 
conditions  be  periodic  and  linear. 

For  constant  thermophysical  properties,  the  principle  of  super¬ 
position  can  be  applied  to  solving  the  heat  conduction  equation. 
Consider  the  following  hypothetical  experiment.  Suppose  that  a  one- 
dimensional  wall  is  initially  at  uniform  temperature  T  ,  and  that  the 
outside  surface  of  the  wall  is  suddenly  raised  to  a  temperature  one 


8 


unit  above  the  initial  temperature,  while  t he  inside  surface  is  main¬ 
tained  at  T  .  The  resulting  heat  flux  at  the  inside  surface  will 

o 

vary  with  time,  as  sfiown  qualitatively  in  Figure  1-1.  Note  that  if 
the  outside  wall  temperature  had  been  raised  by  two  units  instead  of 
one,  the  response  (heat  flux  versus  time)  would  simply  be  twice  that 
shown  in  Figure  1-1.  Similarly,  the  response  to  a  series  of  step 
changes  in  outside  surface  temperature  occurring  at  different  times 
is  obtained  by  adding  the  responses  resulting  from  each  step  change. 
By  representing  the  outside  surface  temperature  as  a  series  of  posi¬ 
tive  and  negative  step  changes  of  appropriate  magnitude  (rectangular 
pulses),  the  flux  caused  by  any  arbitrary  temperature  variation  can 
be  determined  if  the  flux  resulting  from  a  unit  step  change  is  known. 

The  earliest  methods  for  determining  a  multilayered  slab's  heat 
flux  response  to  a  unit  step  temperature  change  were  numerical  (Nessi 
and  Nissole  [ref  12],  Brisken  and  Reque  [ref  13]).  That  is,  for  a 
given  wall  or  roof  section,  the  response  was  determined  once  by  solv¬ 
ing  the  lumped  parameter  resistor-capacitor  analog,  yielding  a  set  of 
so-called  response  factors.  Thereafter,  any  flux  was  determined  by 
applying  the  response  factors  to  the  actual  temperature  profile  as 
approximated  by  rectangular  pulses.  These  methods  were  extended  (by 
analogy  to  electrical  networks)  so  that  the  step  response  of  entire 
rooms  could  be  determined. 


Mitalis  and  Stephenson  (refs  14,15,16),  and  Kusuda  (ref  17)  made 
important  improvements  in  the  response  factor  method  by  approximating 
temperature  as  overlapping  triangular  pulses  (equivalent  to  a  tra¬ 
pezoidal  temperature  profile  approximation)  and  by  treating  the  mul¬ 
tilayered  slab  by  exact  analysis  rather  than  by  lumped  parameter 
methods.  The  approach  consists  of  taking  the  Laplace  transform  of 
the  heat  conduction  equation  and  boundary  conditions  and  solving  the 
equation  in  the  Laplace  domain.  The  required  boundary  condition 
matching  at  the  interface  of  intermediate  layers  of  a  muln layered 
slab  yields  (in  the  Laplace  domain)  an  overall  transfer  function  for 
a  multilayered  slab.  After  multiplying  the  transfer  function  by  the 
Laplace  transform  of  the  triangular  temperature  pulse,  the  inverse 
transform  is  calculated  by  numerically  solving  for  roots  of  the 
characteristic  equation  and  summing  the  residues  at  each  of  these 
poles.  Carrying  out  this  procedure  for  heat  flow  in  both  directions 
yields  three  one-column  matrices  --  [X],  [Y],  and  [Z]  --  called 
response  factors  that  relate  heat  flow  to  surface  temperature  accord¬ 
ing  to  the  following  formula: 


q  .  =  XT  .  •  X •  -  XT...Y. 

°’t  jV-t-J  j  j=o  1  J 


[Eq  1-2] 


q.  =  XT  .  .  Y-  -  XT.  .  •  Z  ■ 
Vt  J=(J  0,t-j  j  j=0i,t-j  j 


[Eq  1-3] 


where  : 


q  .is  the  outer  surface  heat,  flux  at  time  t 
q.  is  the  inner  surface  heat  llux  at  time  t 
Tq  1  .  and  ! .  are  the  current  and  past  outside 
and  inside  surface  temperatures. 

Note  that  as  j  increases,  the  response  factors  converge  to  a 
common  ratio.  This  fact  and  other  algebraic  manipulations  developed 
by  Peavy  (ref  18)  allow  surface  heat,  fluxes  to  be  calculated,  using 
current,  surface  temperatures  and  a  modest  number  (fewer  than  20)  of 
surface  temperature  and  surface  flux  histories.  Chapter  2  provides  a 
thorough  exposition  of  response  factor  methods  for  solving  transient 
multilayered  slab  heat  conduction  problems. 

This  response  factor  conduction  model  is  used  as  part  of  a 
surface-by-surface ,  hour-by-hour  heat  balance  load-predicting  method 
in  the  BLAST  (ref  19)  and  N6SLD  (ref  20)  computer  programs. 

A  final  extension  of  this  method,  also  from  Mitalis  (sec  Appen¬ 
dix  B  of  ref  20),  allows  for  decoupling  of  surface  heal  fluxes  by 
using  weighting  factors.  If  the  heat  transfer  in  a  room  is  assumed 
to  be  linear  (linearized  radiant  heat  transfer)  and  invariant  (con¬ 
stant  tfiermophysu-.il  and  convective  heat  transfer  coefficients),  llun 
f.hf  room  luid  clue  to  a  unit  heat  flux  pulse  can  be  determined  by 
applying  detailed  fie  a  f  balance  computer  techniques  such  as  those  used 
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in  BLAST  or  NBSLIJ.  Weighting  factors  can  then  be  calculated  which 

relate  the  room  load  at  Lime  t,  t+1,  t+2 .  to  a  unit  flux  at  time 

t.  Mitalis  applied  Z-transform  methods  to  the  calculation  of  weight¬ 
ing  factors  for  conduction,  transmitted  solar  radiation,  internal 
radiation,  and  space  temperature  drift  for  light,  medium,  and  heavy 
buildings  (ref  21).  Cumal i ,  et  al .  (ref  22),  have  recently  studied 
the  sensitivity  of  weighting  factors  to  variations  in  room  geometry, 
mass,  and  other  variables. 

Proposed  Method  for  Space  Load  Prediction 

Calculating  hourly  loads  over  a  1-year  period,  although  accu¬ 
rate,  can  be  time-consuming  and  somewhat  costly.  Yet  for  design  of 
building  mechanical  systems,  comparative  analyses  are  usually  more 
important  than  absolute  accuracy.  Hence,  we  will  propose  a  stream¬ 
lined  load-predicting  procedure  which  will  characterize  the  important 
dynamics  of  building  heat  transfer,  but  will  not  require  detailed 
hourly  heat  balance  cal cul ations.  Such  a  procedure  must  be  able  to 
produce  a  time  series  of  predicted  loads  which  capture  the  general 
effects  of  local  climate,  building  use,  and  building  construction 
without  necessarily  predicting  building  response  for  a  specific 
recorded  weather  data  set. 
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1.5 


The  following  combination  and  extension  of  earlier  harmonic  and 
weighting  factor  methods  appears  to  minimize  required  calculations, 
yet  still  considers  dynamic  response. 

1 


Using  hourly  cl imatol og ical  records  (on  computer  tape)  and  fast 
Fourier  transforms  or  other  time  series  analysis  techniques,  deter¬ 
mine  the  most  significant  coefficients  and  frequencies  of  a  Fourier 
series  approximation  to  weather  data  variables.  A  set  of  coeffi¬ 
cients  will  be  required  for:  (1)  dry-bulb  temperature,  (2)  sol-air 
temperature  for  surfaces  of  various  solar  absorptivities  and  facing 
directions,  (3)  incident  solar  radiation  on  windows  facing  different 
directions,  and  (4)  dew  point  or  specific  humidity. 

Step  P. 

For  each  different  wall  and  roof  construction  and  for  each  dif¬ 
ferent  exposure  and  for  each  different  frequency  judged  to  be  impor¬ 
tant  in  step  1,  determine  the  conductive  heat  flux  resulting  from 
unit  sinusoidal  and  cosinusoidal  sol-air  temperature  variation. 

Step  3 

For  each  exposure  and  frequency,  multiply  the  conductive  heat 
flux  coefficients  from  step  2  by  the  area  of  the  surface  and  then  by 
the  sol-air  temperature  coefficients  determined  from  step  1.  This 
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will  provide  the  frequency  domain  representation  of  conductive  heat 
flux  into  the  space  (assuming  the  space  is  at  constant  temperature). 

Step  •! 


Evaluate  the  function  obtained  in  step  3  at  equally  spaced 
points  in  time.  In  addition,  determine  the  solar  heat  gain  transmit¬ 
ted  through  windows  and  the  scheduled  heat  gain  from  people,  lights, 
and  equipment. 

Step  5 

Apply  appropriate  weighting  factors  to  each  heat  gain  and  sum  to 
determine  space  loads. 

The  following  very  simple  example  illustrates  this  method. 

Suppose  a  room  lias  a  south-facing  wall  of  area  and  a  west- 
facing  wall  of  area  A  The  rest  of  the  room  surfaces  are  interior 
partitions,  a  ceiling,  and  a  floor.  (Assume  a  second  story  above  the 
room  and  a  basement  below.)  Suppose  that  from  step  1  of  the  above 
procedure,  the  sol-air  temperature  for  south-facing  surfaces  with 
absorptivity  corresponding  to  the  absorptivity  of  the  south  wall  of 
the  room  can  be  approximated  by: 

sol-air  south  Oss  Is  '  1  '  Is  '  1 


+  S365sSin^36b^  +  C36bs  COS^w3bSt^  ^Eq 


where  t  is  time,  and  the  Cs  and  Ss  are  constant  coefficients 


1  0 


corresponding  to  important  weather  frequencies,  w  (the  mean,  one 
cycle  per  year,  and  365  cycles  per  year  in  the  above  example).  Simi¬ 
larly,  suppose  that  the  sol-air  temperature  for  the  west  wall  is: 


Tsol-a,r  west  ‘  C0w  *  Slwsin<“lt>  *  tl,cos<ult) 

*  s365wsin(aW>  *  L365wcos(“36bti  ^ 

If,  from  step  2,  the  heat  flux  coefficients  for  the  two  walls  at  the 
vanous  important  frequencies  are  UQ,  U1$s,  URs,  U365ss>  U36bcs, 

U0w’  UIsw’  Ulcw’  U365sw’  and  U365cw’  then  the  heat  gain  through  thG 
exterior  walls  as  a  function  of  time  is: 

Vl>  ■  \U0s  <C0s  -  V  *  Vow  <C0w  '  TR> 

*  <\  U1SS  Sls  *  A,  Ulsw  Slw)Si"<“ll> 

*  (AS  “lCS  C1S  *  AW  UlCW  Clw)C0S(U'l!') 

*  <As  U365ss  S365s  4  Aw  U365ws  WS1 n*“366l * 

*  <As  U365cs  C 365s  *  ftw  U365cw  C365w>cos(“365t'  [Eq 

This  function  can  be  evaluated  at  any  time  t  to  determine  the  heat 
flux  through  exterior  walls  (step  4). 
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The  heating  or  cooling  load  in  the  room  at  time  t  for  the 
assumed  constant  room  temperature,  Tp,  is  obtained  by  adding  the  con¬ 
tributions  of  wall  conductive  flux,  solar  energy  transmitted  through 
glass,  lighting  energy,  and  the  energy  emanating  from  people,  with 
each  weighted  according  to  step  5.  For  example,  if  the  lighting 
energy  flux  into  the  room  is  (t),  then  the  contribution  to  room 
load  from  the  lights  is: 

LL(t)  =  V0  QL(t)  +  V1  QL(t_At)  +  V2  QL(t-2At)-.. 

-  W}  LL(t-At)  -  W 2  LL(t-2At)  -  W3  LL(t-3At)...  [Eq  1-7 j 

where  Q^(t-  n  At)  are  past  values  of  lighting  flux,  and  L^(t-  n  At) 
are  past  contributions  of  lighting  energy  to  the  space  load. 

The  total  load  at  time  t  is  therefore: 

m 

L(t)  =  l  U(t)  [Eq  1-8] 

k=l  K 

where  the  sum  is  taken  over  the  number  of  heat  gain  components. 

The  loads  from  step  5,  along  with  a  matching  weather  data  time 
series  (obtained  by  inverting  the  transforms  obtained  in  step  1), 
will  provide  the  input  necessary  to  simulate  air  distribution  systems 
and  boiler/chiller  plants. 
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The  proposed  load  calculation  procedure  has  several  advantages 
over  existing  hour-by-hour  methods.  The  number  of  calculations 
required  is  reduced  substantially,  particularly  since  step  1  --  the 
processing  of  weather  data  —  need  only  be  accomplished  once  for  any 
given  location.  The  major  part  of  the  load  calculation  procedure  -- 
the  evaluation  of  sine  and  cosine  terms  --  is  simple  enough  to  be 
implemented  on  modern  microcomputers  (for  example,  as  a  subprogram  it 
a  larger  system  simulation  program).  This  should  substantially 
enlarge  the  number  of  engineers  and  architects  having  the  resources 
and  skill  to  apply  the  method  routinely.  The  graphic  portrayal  of 
the  frequency  response  of  wall  and  root  sections  may,  in  itself, 
prove  to  be  a  valuable  design  aid. 

The  proposed  method  also  has  limitations.  First,  the  approxima¬ 
tion  of  weather  data  by  a  few  terms  of  a  Fourier  series  necessarily 
"filters  out"  climate  extremes.  Flence,  energy  use  estimates  based  on 
this  filtered  data  can  be  expected  to  be  low  (ret  2).  Second,  the 
methou  is  less,  rather  than  more,  detailed  and  therefore  suffers  from 
at  least  as  many  simplifying  assumptions  as  the  hour-by-hour  method. 

Object ive 

The  objective  of  this  research  was  to  develop  the  necessary 
mathematical  techniques  to  implement  and  test  the  new  method  for 
space  load  prediction  described  in  the  previous  section.  Because 
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much  of  the  complexity  of  space  load  calculation  stems  from  the  need 
to  deal  with  transient  heat  flow  through  multilayered  slabs,  a  second 
objective  was  to  provide  a  comprehensive  mathematical  development  of 
response  factor  theory. 

Approach 

Four  main  tasks  were  performed  to  achieve  the  objective: 

1.  A  procedure  was  developed  to  calculate  the  frequency 
response  of  multilayered  slabs. 

2.  Various  procedures  were  applied  to  the  analysis  of  weather 
data,  and  the  best  method  for  determining  the  periodic  components  of 
the  data  was  selected. 

3.  Errors  associated  with  the  use  of  synthesized  weather  data 
which  contained  only  sinusoidal  components  were  assessed. 

4.  The  proposed  method  for  space  load  prediction  was  tested. 


Outl ine 


The  results  of  this  research  are  described  sequentially  in 
Chapters  2  through  7.  Chapter  2  provides  a  detailed  development  of 
transfer  function  methods  applied  to  heat  conduction  in  multilayered 
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slabs.  It  provides  the  necessary  background  for  the  study  of  the 
frequency  response  of  rnul ti 1 ayered  slabs  described  in  Chapter  3. 

Chapter  4  describes  a  thorough  analysis  of  important  weather 
data  variables.  Chapter  5  describes  tests  performed  to  determine  the 
error  in  estimated,  annual  building  energy  use  associated  with  using 
only  the  sinusoidal  components  of  weather  data  instead  of  actual 
weather  data.  Chapter  6  presents  an  hour-by-hour  comparison  between 
results  of  the  proposed  load  calculation  method  and  more  detailed 
methods. 

Chapter  7  summarizes  the  results  of  the  research  and  outlines  a 
procedure  for  implementing  the  new  load-predicting  procedure  in  a 
microprocessor-based  building  energy  analysis  program. 


1 

t 
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2  ONI -DIMENSIONAL  HEAT  CONDUCTION  THROUGH  MULTILAYERED  SLABS 


The  Heat  Conduction  Equation 


Heat  conduction  through  a  one-dimensional  homogeneous  slab  is 
governed  by  the  following  second-order  partial  differential  equation 
(see  Figure  2-1 ) : 


a2T(x,t)  =  1  3T(x,t) 
i)x2  01  3t 


[Eq  2-1] 


where  T  is  the  temperature  at  position  x  and  time  t,  and  <*  is  the 

i 

thermal  diffusivity  a  =  — where  k  is  the  thermal  conductivity 

(W/m°K),  Pis  density  (kg/m3),  and  Cp  is  the  specific  heat  (J/kg°K). 
The  heat  flux  at  any  position  x  and  time  t  is  given  by: 

q(x,t)  =  -k  .  [Eq  2-2] 

In  both  the  above  relations,  k,  p,  and  Cp  were  assumed  to  be  con¬ 
stant  . 


A  common  approach  to  the  solution  of  the  above  equations  is  to 
use  the  Laplace  transform,  which  is  defined  for  any  transformable 
function  f(t)  as 


A[f(t)]  =  F(s)  = 


I 


f(t)  e"stdt 


[Eq  2-3] 


0 


where  s  is  a  complex  number.  The  utility  of  the  Laplace  transform 
stems  from  the  following  property: 


C^]  =  sF(s)  -  f(t) 


t=0 


[Eq  2-4] 


Hence,  the  Laplace  transform  of  Equation  2-1  transforms  this  partial 
differential  equation  into  an  ordinary  differential  equation: 


dJLIx^l  -  I  sT(x>s, 


dx 

where  we  have  assumed  that  T(x,0)  =  0.  The  solution  of  this 
transformed  differential  equation  is 


[Eq  2-5] 


T(x,s)  =  A  cosh(x  /s/a)  +  B  sinh(x  /s/a) 


[Eq  2-6] 


The  transform  of  Equation  2-2  is: 

q(x>s)  =  -k  [Eq  2-7] 

where  we  have  again  assumed  that  T(x,0)  =  0.  On  differentiation  of 
Equation  2-6  with  respect  to  x  and  substitution  into  Equation  2-7,  we 
have : 


q(x,s)  =  -k  /s/a  A  sinh(x  /s/a)-  k  /s/a  B  cosh(x  /s/a)  [Eq  2-8] 
If  we  now  consider  only  ttie  temperature  and  heat  flux  at  the 


surfaces  of  the  slab  (at  x  =  0  and  x  -  1  where  l  is  the  thickness  of 
the  slab  in  meters),  we  can  write  from  Equations  2-6  and  2-8: 


T-j  (s  )  =  A  [Eq  2-9] 

q-j  ( s )  =  -kB  /s/a  [Eq  2-10] 

T,,(s)  =  A  cosh(i  /s/u)+  B  sinhU  /sf  •. )  [Eq  2-11] 

q2(s)  =  -k  /s/a  A  sinh(£  /s/a)-  k  2? Ju  B  coshU  /sTa)  [Eq  2-12] 


where  T^(s)  =  T(0,s),  which  is  the  transform  of  the  temperatures  at 
surface  1 

( s )  =  q(0,s),  which  is  the  transform  of  the  flux  at 
surface  2 

T^s)  and  q2 ( s )  are  corresponding  transforms  of  temperature 
and  flux  at  surface  2  [i.e.,  T(£,  s)  and  q(£,  s)]. 

We  may  now  specify  any  two  of  T^(s),  q^(s),  T^s),  or  q^s)  as  boun¬ 
dary  conditions  (transforms  of  the  appropriate  time-dependent  boun¬ 
dary  conductions)  in  order  to  complete  the  solution.  While  the  phy¬ 
sical  problem  will  dictate  which  boundary  conditions  are  known,  we 
will  temporarily  assume  that  both  T^s)  and  q2(s)  are  known. 
Proceeding  with  this  assumption,  we  can  use  Equations  2-11  and  2-12 
to  find  A  and  B  in  terms  of  the  physical  constants  of  the  slab  and 
the  Laplace  transform  parameter,  s. 


[Eq  2-1 3] 


( s )  -  B  sinb(?  /s/a) 
COSh(i  /s/a) 


-qp(s)  -  k  /s/a  A  sinh  </,(/ s/a) 

B  =  - - . . -  [Eq  2-14] 

k  /s/u  cosh(k  /s/a) 

A  few  algebraic  manipulations  and  hyperbolic  trigonometric  identities 
yi  cl  d : 

A  =  (cosh(£  /i7a))T„(s)  +  /S1nh(^  A/a)\  q  (s)  [Eq  2-15] 

*  V  k  /s/^  /  c 


B  =  (-sinh(£  /s/a) )  T?(s)  -  -s-/ct-i-  [Eq  2-16] 

k  /s/a 

Substituting  for  A  and  B  in  Equations  2-9  and  2-10  yields  the 
required  solution  in  the  Laplace  domain: 


T,(s)  =  (cosh(£  /s/a))  T,(s)+  ( — ^=-sinh(£  / s/a))  q~(s )  [Eq  2-17] 

1  1  \k/s/a  /  * 


q  1  ( s  )  =  (k/s/a  sinh(£  /s/a^Tgls)  +  (cosh(£  /s/a))q2(s)  [Eq  2-18] 


For  notational  convenience,  we  may  now  define  new  variables  as 


fol 1 ows : 


f 


A  ( s ) 

=  cosh  ( v  /s T(t) 

B(s ) 

=  — - --  sinh  ( '.■ 

VS./  t) 

k/s/u 

C(s) 

=  k  /s  /  a  s  i  n  h  ( , 

2s/0 

D(s) 

=  cosh ( ,  Af  t) 

where  A(s)  and  B(s)  should  not  be  confused  with  A  and  B  used  previ¬ 
ously.  With  these  new  variables.  Equations  2-17  and  2-18  become: 


T-,  ( s )  =  A(s)  T2(s)  +  B(s )  q9(s) 


q1  (s)  =  C(s)  T2(s)  +  D(s)  q2(s) 


[Eq  2-19] 

[Eq  2-20] 


Note  that  these  two  equations  can  be  solved  for  any  two  unknowns  in 
terms  of  the  two  knowns  which  arise  from  the  physical  boundary  condi¬ 
tions  of  the  problem.  For  example,  if  the  transforms  of  temperature 
variation  with  time  are  known  for  both  surfaces,  then  the  fluxes  are: 


(s) 


q2(s) 


Tl(s) 


BTsT  T1(S) 


eT[s7  Vs)’ 


B(s )  t  0 
[fj  T2(s),  B(s )  t  0 


[Eq  2-21] 


[Eq  2-22] 


At  least  in  theory,  if  Laplace  transforms  exist  for  the  two 
required  boundary  conditions,  the  unknown  fluxes  or  temperatures  can 
be  determined  by  finding  the  inverse  Laplace  transforms  of  equations 


like  2-19,  2-20,  2-21,  or  2-22. 

We  may  now  extend  the  treatment  of  one-dimensional  heat  flow  to 
include  multilayered  slabs.  Notice  that  Equations  2-19  and  2-20 
describe  the  transform  of  heat  flow  and  temperature  on  one  surface  in 
terms  of  the  transform  of  the  heat  flow  and  temperature  on  the  other 
surface.  We  can  rewrite  these  equations  in  matrix  form  as: 


Vs  >1 

~A(s) 

B(s)“ 

"t2(sT 

9^ (s )_ 

_C(s) 

D(s)_ 

_  92  ( s )  _ 

[Eq  2-23] 


Suppose  now  that  we  have  a  two-layer  slab  (see  Figure  2-2),  where 
Tj(s)  and  q^(s)  are  the  transforms  of  temperature  and  flux  at  the 
inner  surface  of  slab  1,  and  T^s)  and  q^s)  refer  to  the  interface 
between  slabs  1  and  2.  T^(s)  and  q^s)  refer  to  the  outer  surface  -- 

surface  3.  We  may  treat  each  slab  individually  and  apply  Equation 
2-23  to  each.  Noting  that  surface  2  is  both  the  "outside"  of  slab  1 
and  the  "inside"  of  slab  2,  we  have  for  slab  1: 


SLA8  ( 

SLAB  2 

T' ! 

M 

M 

1 

J 

U  ) 

x  W 

*  I 

f  ~ 

j 

1 

SURFACE  / 

•  l\  ^ 

•  iZ  m 

^  SURFACE  Z 

^  SURFACE  3 

Figure  2-2.  One-dimensional  two-layer  slab 


[Eq  2-24] 


Vs)- 

_A-j  ( s ) 

B-j  ( s  r 

'T2(s)~ 

_q  i  ( s )_ 

_C-j  (s) 

D]  (5  )_ 

_q2(s)_ 

and  for  slab  2: 


“t2(s)" 

_A2(s) 

B2(s) 

"T3(s)- 

_q  2  ( s )_ 

_C2(s) 

D2(s)_ 

_q3(s)_ 

where 

A-j(s)  =  coshU-j  /s/a-j ) 

B-j  ( s )  =  (l/k^»^s7a^)  sinh(v.1  /s/o^ ) 
Cj(s)  =  k.j  /s/«i  (sinh(£-|  /s/a^) ) 
(s  )  =  COShUi  /s/a-|  ) 


are  all  dependent  on  the  properties  of  slab  1  and  the  Laplace 
transform  variable,  s.  A2(s),  B2(s),  C2(s),  and  D2 ( s)  are  similarly 
defined  based  on  the  properties  of  slab  2. 


We  may  now  simply  substitute  the  right-hand  side  of  Equa- 

FM 


tion  2-25  for 


q2(s) 


in  Equation  2-24,  yielding 
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'Vsr 

A-j  (s) 

B-j  (s)~ 

-qi  (sL 

_C i  ( s ) 

D-j  ( s )_ 

ft2(s)  &2(s) 


_C2(s)  D?  ( s )_ 


-T3(s)- 

_q3(s)_ 

[Eq  2-26] 


The  recursive  extension  of  Equation  2-26  to  an  n-layered  slab  is  now 
obvious : 


"T^sf 

p 

_Vs)_ 

l< 

A-j  (s)  B-j  (s ) 


A3(s)  B3(s) 


Lq1(s)J  1C]  (s)  D^sy  ic2(s)  D2(s)J  [c3(s)  D3(s)J 


An(s)  Bn(s)' 


C  (S)  D  (s) 


"a2(s)  e 

_ 

_C2(s)  I 

- 

~Ws)~ 

-Vl (s)- 

[Eq  2-27] 


This  fundamental  equation  permits  the  calculations  of  the  so-called 
"transmission  matrix,"  defined  as 


3' 


A(s) 

B(sf 

\  (s) 

B  -|  ( s  )_ 

A2(s) 

B2(s)' 

C(s) 

D(s)_ 

C-j  (s) 

D1  ( s )_ 

_C2(s ) 

D2(s)_ 

[Eq  2-28] 


Vi(s>  W5)1 

A„(s)  Bn(s)' 

-Cn<s> 

With  A(s) ,  B(s) ,  C(s) ,  and  D(s)  now  redefined  as  elements  of  the 
above  overall  transmission  matrix,  the  multilayered  slab  problem 
takes  the  same  form  as  that  of  a  single-layer  slab.  Equations  2-19, 
2-20,  2-21,  and  2-22,  for  example,  can  be  written  for  a  multilayered 
si ab  as : 


T i  ( s  )~ 

~A(s) 

B(s)~ 

Vi<s>' 

(s)_ 

C(s) 

D(s)_ 

[Eq  2-29] 
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Ms)  ~ 

IWS> 

D(s)  1 
Bfs)  '  BTi7 

1  A(s) 

erry  "  b(sT- 


"Vs) 

-Ws>- 

[Eq  2-30] 


While  Equations  2-29  and  2-30  are  conceptually  simple,  their  practi¬ 
cal  use  is  mathematically  tedious  for  all  but  the  simplest  of  cases. 
We  will  therefore  first  consider  single-layer  slab  applications  and 
then  extend  these  results  to  multilayered  slabs. 

Before  proceeding  with  an  example,  let  us  note  some  of  the  pro¬ 
perties  of  the  transmission  matrix  for  a  single-layer  slab. 


"Vs) 

8,(5)" 

rcosh(£.  /s/a,),  sinhU.  /s/a,) 

k,/s/a,  ‘ 

Vs) 

D,(s)_ 

\k]  /s/a,  sinhU,  /s/a,),  coshU,  /sTa^lJ 

[Eq  2-31] 


We  can  also  write  this  matrix  in  terms  of  the  thermal  resistance,  Rj , 
and  thermal  capacitance,  C, ,  defined  as  follows: 


R 


1 


and,  in  general ,  Rn 


[Eq  2-32] 


C1  =  VlCPl  and’  general,  Cp  =  £nPnCpn 


[Eq  2-33] 
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£  n 

Notice  that  R  C  =  — .  Thus,  the  single-layer  transmission  matrix 
n  Pn  01  n 

becomes : 


"A1  (s) 

B  i  ( s  )~ 

_C  t  ( s ) 

D1  (s)_ 

cosh  /s R-j  C'i"  , 


/sR1C] 


/s  R-j  C-| 


sinh  /sR-jC^j" 


[Eq  2-34] 


(_ — ^ — L  sinh  /sR-jC-j  ,  cosh  /sR^C^J 


The  determinant  of  this  matrix  is  cosh  (/sR^Cj )- sinh  (/sR^C^)  =  1. 

The  multilayered  slab  transmission  matrix  therefore  also  has  a  deter¬ 
minant  of  one.  Also,  as  the  capacitance  of  the  material  goes  to 


ySRici  J  Rj 

zero,  cosh/sR^  =  1  and  — p — -si nhi/sR^Cj"  =  0.  sp-  c~- ■ 

can  be  expanded  to  reveal  its  limit  as  -*■  0. 


sinh  / sR^cY 


✓sR]C1 


sinh  /sR-jC-j  = 


/sR1C] 


/srTcT  /srTc, 
/sR,C.  +  - —  +  - — 


_5 


'ri 


3! 


5! 


[Eq  2- 


R,  _ 

As  C.  -Ml,  — -  "inhv^sRTC,  =  R 

1  /sTjiq  1  1 


A-j  (s )  B-j  (s ) 

1 

Rl‘ 

.C^s)  D-j  ( s )_ 

_0 

1  _ 

Thus,  for  a  1 i ghtwei ght  slab: 


[Eq  2-36] 


This  means  that  air  layers  or  other  light  material  can  be  routinely 
included  in  the  calculation  of  the  multilayered  transmission  matrix. 
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Lot  us  now  consider  the  case  of  a  single  layer  exposed  on  one 
side  to  a  temperature  "ramp"  while  the  temperature  on  the  other  side 
is  fixed.  This  case  is  of  practical  interest,  because  an  arbitrary 
surface  temperature  variation  can  be  approximated  by  a  series  of  such 
ramps  with  alternating  positive  and  negative  slopes.  Figures  2-3  and 
2-4  show  how  three  ramps  can  be  used  to  form  a  triangular  pulse  and 
how  overlapping  triangular  pulses  form  a  trapezoidal  approximation  to 
an  arbitrary  surface  temperature. 


The  Laplace  transform  of  a  ramp  with  unit  slope 

if  this  transformed  boundary  condition  is  applied  to 
single-layer  slab  while  the  temperature  on  side  1  is 
Equation  2-25  (or  Equation  2-30)  becomes: 


is  —.  Hence, 
s^ 

side  2  of  a 
held  at  zero, 


q]  (s) 

“D1  (s) 

1 

0 

‘  b^TsT 

q2(s) 

i 

Vs) 

l 

Lb^T 

VsT_ 

Ls2  J 

[Eq  2-37] 


If  we  are  interested  in  the  conductive  flux  at  surface  1  due  to  the 
ramp  boundary  condition  on  surface  2,  we  have 


q^s)  =  - 


/s  R-|  C-j 

Rl  sinh  /sRi  C-j 


[Eq  2-38] 


TEMPERATURE 


Figure  2-4.  Overlapping  triangular  pulses  as  an  approximation 
to  a  continuous  function. 


1 


The  inverse  Laplace  transform  yields  flux  variation  at  surface  1  as  a 
function  of  time. 


1  [q-j  (s )  ]  =  q-j(t)  =  :  1 


/sR^" 


R-|  sinh 


[Eq  2-39] 


Applying  the  general  formula  for  finding  inverse  Laplace 
transforms,  we  have: 


q-|  ( t ) 


1 

2irj 


a+jaj 

q-|  (s)eStds 

a-j® 


[Eq  2-40] 


where  j  =  /T  and  a  is  a  very  large  real  number.  From  complex  vari¬ 
able  theory,  the  above  integral  is  equal  to  the  sum  of  the  residues 
at  the  poles  of  qj(s)est.  Poles  are  points  where  s  assumes  a  value 
so  as  to  make  q^(s)es^  undefined.  In  our  case,  any  time  the  denomi¬ 
nator  of  the  right  side  of  Equation  2-38  is  zero,  the  function  is 

undefined.  Hence,  poles  exist  at  s  =  0  and  at  /sR^C^  =  nnj  or 
2  2 

s  =  -n  tt  /RjCj,  where  n  =  1,2,3,....  Notice  that  the  second-order 

pole  at  zero  is  due  to  -L,  the  boundary  condition,  and  that 

s^ 

sinh  / sR^Cj  does  not  represent  a  pole  at  zero,  since,  as  we  have 

/sTTcT 

already  shown,  lim  =  1.  All  other  poles  lie  on  the 

s  -  0  sinhvSR1C1 

negative  real  axis  in  the  s-plane. 


which  is  the  conductive  heat  flux  at  surface  1  at  time  t  due  to  a 


ramp  temperature  increase  of  unit  slope  applied  to  surface  2  at  time- 
zero  while  maintaining  surface  1  at  zero  temperature. 


We  are  now  ready  to  calculate  the  flux  due  to  the  triangular 
pulse  of  Figure  2-3.  The  flux  at  time  5  is: 


,  [C,  *  ~  (-1  )nC1e'Sn2l,2/RlCl 

%  6  -  nTT 


[Eq  2-45] 


where  we  have  simply  replaced  t  by  6  in  Equation  2-41  and  replaced 
the  unit  ramp  with  a  ramp  of  slope  T/6.  At  t  =  2  6  ,  both  ramp  1  and 
ramp  2  of  Figure  2-3  contribute  to  the  flux.  The  contribution  of 


ramp  1  is: 


q, (25)1 


due  to 
ramp  1 


f>  ,  nnr  - 26 n^Tr  / R, C, 

_  T  C1  26  ,  0  "  ('1}  Cle  1  1 

6  6  R.  d  4  2  2 

1  n=l  n  tt 


[Eq  2-46] 


Ramp  2  starts  at  time  6  and  has  a  slope  of  -2T/6.  Hence,  its  contri¬ 
bution  to  the  flux  is: 


(2$) 


due  to 
ramp  2 


fc  oo  (.1  inr  e-(26-6)n  it  /R,C. 

2T  h  26-6  A  0  ”  C1  1 

6  6  “  R,  S  2  2 

I  n=l  n  it 


[Eq  2-47] 


The  principles  of  superposition  allow  us  to  add  these  two  contribu¬ 
tions  to  yield 


q-,  (26)  = 


+  2 


c 

L 

n  =  l 


(-1 )nc. 


.  2e 


2  2 
-<in  n 


-  e 


2  2 


25n2ir2/R1C1j 


[Eq  2-46] 


Applying  the  same1  p -ocedures  for  t  =  3 6  and  including  t fie  effect  oi 
ramp  3,  which  starts  at  25  and  has  a  slope  of  T/5,  we  have: 


q-,  ( 35 ) 


2  l 

n=l 


(-lA, 


-36n2iT2/  R-j  C-j 


?  2  2  2 
-24n"ir  /R]C1  -6n  V /RiCj 

2e  +  e 


~TT 

n  it 


')] 

I 

J 


[Eq  2-49] 


Similarly,  for  any  incremental  time,  m  6  ,  where  m  is  greater  than  or 
equal  to  3: 


q-|  (nv5 )  = 


T 
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-(m-1 )56  - (m-2)5B  • 

2e  n  +  e  n 


[Eq  2-50] 


where  B 


2  2/d 
n  it  /R 
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Equations  2-45,  2-48,  and  2-50  can  be  used  together  to  calculate 
the  flux  at  surface  1.  Figure  2-5  illustrates  typical  variations  in 


flux  with  time  for  a  heavy  and  light  slab  subjected  to  a  triangular 
pulse  starting  at  t  =  0.  Note  that  flux  is  defined  to  be  positive  if 


it  flows  in  the  positive  x  direction;  hence,  q^  is  the  heat  flux  out 
of  surface  1. 


Suppose  now  that  we  wish  to  calculate  the  conductive  flux  at 
surface  1  due  to  a  change  in  the  temperature  of  surface  1  itself, 
while  surface  2  is  held  at  constant  temperature.  Following  the  same 
procedure  of  approximating  the  temperature  as  a  series  of  ramps. 
Equation  2-30  becomes: 


[Eq  2-51] 


Substituting  for  D^(s)  and  B^(s)  yields  the  counterpart  to  Equation 
2-38. 


q1  (s) 


i 


1  (cosh  /sR-|C.|)  /sR^C^ 
s^  R-j  sinh  /sR^C^ 


[Eq  2-52] 


2  2 

The  poles  of  this  function  are  again  at  s  =  0  and  s  =  -n  tt  /R^C^. 
Applying  the  inverse  Laplace  transform  formula  and  the  residue 
theorem  as  before,  we  have: 


1  i/i) 
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[Eq  2-53] 
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[Eq  2-54] 


The  counterpart  to  Equation  2-44  is  therefore 


2  2 

-tn Sr  /R,C, 

C1  t  °°  Cle 

ql (t)  =  T  +  R7  '  2  ^  2~2 

!  n=l  n  tt 


[Eq  2-55] 


which  is  the  conductive  heat  flux  at  surface  1  at  time  t  due  to  a 
ramp  temperature  increase  of  unit  slope  applied  to  surface  1  at  time 
zero  while  maintaining  surface  2  at  zero  temperature. 


The  flux  due  to  a  triangular  pulse  applied  to  surface  1  is 
obtained  in  the  same  manner  as  the  derivation  of  Equations  2-45 
through  2-50.  Hence: 
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[Eq  2-56] 
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-26n2TT2/R,C1  -6n2iT2/R,C 
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r  r  f  'm6B" 

T  00  Ci  ie 

a^nv1,)  =  -  g  2  £ - 

1  6  [_  n=l 

2  2 

for  m  _>  3  and  =  n  tt  /  R-j  C-j  . 

Applying  these  equations  produces  Figure  2-6,  the  counterpart  to  Fig¬ 
ure  2-5,  showing  variations  in  conductive  flux  surface  at  1  when  a 
triangular  temperature  pulse  is  applied  to  surface  1  (surface  2  is  at 
zero  temperature). 


Response  factors 


Let  us  abbreviate  Equations  2-56,  2-57,  and  2-58  by  defining  the 


following  variables:* 


»  .i  i  +  i  ?  v  ii! 

A1  "  6  3  R  "  c  S  2  2 
1  n=l  n  it 


[Eq  2-59] 
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[Eq  2-60] 
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[Eq  2-61] 
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,  m  =  3,4,5 


Next  we  abbreviate  Equations  2-45,  2-48,  and  2-50  by  defining  the 
following  variables. 


*  In  some  of  the  literature,  X.,  X^,...  may  be  defined  starting  with 
the  subscript  zero,  i.e.,  Xg,  Xj,  X^... 
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[Eq  2-62] 
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[Eq  2-63] 
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[Eq  2-64] 
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,  m  =  3,4,5 


Notice  that  we  have  introduced  a  sign  change  in  defining  Yj,  and 


We  may  now  apply  the  superposition  principles  to  find  the  con¬ 
ductive  flux  at  surface  1  due  to  the  combined  effect  of  a  triangular 
pulse  of  height  T.  applied  to  the  inside  surface  (surface  1)  and  a 
triangular  pulse  of  height.  Tq  applied  at  the  outside  surface  (surface 
2),  where  both  pulses  start  at  time  zero  (their  apex  is  at  t  =6). 
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q^fi)  =  Tix1  -  V]  [Eq  2-65] 

q-j  (26)  =  TiX2  -  ToX2  [Eq  2-66] 

q7 (m6)  =  T.Xm  -  ToYm,  m  =  3,4,5 .  [Eq  2-67] 


Remembering  that  q^  is  the  flux  in  the  positive  direction,  the  above 
equations  define  the  flux  1 eaving  surface  1  due  to  conduction. 

Ey  simply  changing  our  point  of  reference,  we  can  write  an  equa¬ 
tion  for  the  flux  leaving  surface  1  in  terms  of  present  and  past  sur¬ 
face  temperature  pulses.  Suppose  we  are  sitting  at  a  point  in  time 
coincident  with  the  apex  of  inside  and  outside  triangular  temperature 
pulses  of  height  T-  t  and  TQ  t,  respectively.  Part  of  the  flux  leav¬ 
ing  surface  1  is  due  to  the  upward  part  of  the  ramps  with  slopes 
T.  /6  and  T  ,/6.  This  contribution  is: 

I  *  L  v  j  l 


Vl>t  •  Ti,txl  •  To,tvl 


[Eq  2-68] 


Another  contribution  is  due  to  two  of  the  three  ramps  which  make  up 
the  pulses  centered  at  t-1  (time  is  in  units  of  6  ) .  This  contribu¬ 
tion  is: 


ql ( t)f-l 


Ti ,t-lX2 


T*  Y 
o,t-l  2 


[Eq  2-69] 


Similarly,  all  three  ramps  of  the  pulses  centered  at  t-2  contribute 
to  the  flux  at  t  as  do  the  ramps  making  up  the  pulses  at  t-3,  t-4, 
t-5,  and  so  on.  The  sum  of  all  these  contributions  gives  the  flux 
leaving  surface  1: 


i  ,t 


qj(t)  = 


m=l 


Ti ,t-m+lXm 
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[Eq  2-70] 


The  infinite  series  of  Xs  and  Ys  above  forms  two  parts  of  a  set 
of  three  series  known  as  X,Y,Z  response  factors.  Although  we  have 
only  shown  how  X  and  Y  response  factors  are  calculated  for  a  single- 
layer  slab,  they  are  equally  well  defined  for  multilayered  slabs. 
Repeating  Equation  2-30  will  help  make  their  definition  clear. 


”q1  (s) 
-Vi(s>- 
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- 

[Eq  2-71] 
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The  X  series  of  response  factors  is  defined  as  the  inverse  Laplace 
transform  of  the  quantity  D(s)/B(s)  times  the  Laplace  transform  of  a 
triangular  pulse  of  unit  height.  We  have  al ready  shown  that  this 
pulse  is  made  up  of  three  ramps,  which  leads  to  a  three-part  form  for 
'he  nquirtd  inverse  transform.  If,  for  notational  convenience,  we 
d( tint  the  transform  of  the  pulse  as  P(s),*  we  can  formally  define 
*  f  t  set  of  X  response  factors  as: 


t=m6,  m=l ,2,3. . . . 


[Eq  2-72] 


tins  series  is  sometimes  referred  to  as  the  set  of  internal  response 
factors,  since  it  represents  the  flux  response  of  the  inside  surface 
to  changes  in  inside  surface  temperature. 


The  Y  series  of  response  factors  is  defined  as 


*  The  Laplace  transform  of  a  triangular  pulse  of  unit  height  and 
base  2  6  is: 

P(s)  =  — ~  0  <_  t  ^  6 

6s 

P(s)  =  (l-2e's6)/6s2  6  <  t  <  26 


P(s)  =  (l-2e's*  +  e'2s,')/6s2  t  >  26 


We  could  equally  well  have  found  the  inverse  of 
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[Eq  2-73] 


This  series  is  sometimes  referred  to  as  the  set  of  cross  response 
factors,  since  it  characteri zes  both  the  flux  response  of  the  inside 
surface  to  outside  temperature  variation  and  the  flux  response  of  the 
outside  surface  to  inside  temperature  variation. 


Finally,  the  Z  series  of  response  factors  is  defined  as: 


Z 


m 


t=m6,  m=l ,2,3. . . . 


[Eq  2-74] 


This  series  is  sometimes  referred  to  as  the  set  of  external  response 
factors,  since  it  represents  the  flux  response  of  the  outside  surface 
to  changes  in  the  outside  surface  temperature.  For  a  single-layer 
slab  or  a  symmetric  multilayered  slab,  A(s)  =  D(s)  and  therefore,  Z[n 
=  Xm<  In  all  other  cases  they  are  not  equal. 

With  this  definition  of  the  Z  response  factor  series,  we  can 
develop  an  equation  for  qn+^  (t),  the  outside  flux,  in  exactly  the 
same  way  we  derived  Equation  2-70.  The  final  result  is: 

oo  00 

qo,t  '  qn+l  ^  '  J:-|Ti  ,t-m+lYm  "  To,t-m+lZm 


[Eq  2-75] 


This  is  the  positive  X-directed  flux  at  the  outside  surface.  It  is 
therefore  the  conductive  flux  into  the  outside  surface. 

We  now  take  note  of  the  complexity  involved  in  finding  response 
factors  for  a  multilayered  slab.  Equations  2-72,  2-73,  and  2-74  have 
a  deceivingly  simple  form.  However,  for  a  multilayered  slab,  each  of 
the  terms  A(s),  B(s),  C(s),  and  0 ( s )  represents  complicated  sums  and 
products  of  hyperbolic  functions  of  s  and  the  properties  of  each 
layer.  Recall  particularly  that  we  must  find  the  poles  of  the  vari¬ 
ous  Laplace  transformed  expressions  in  order  to  find  the  inverse 
transforms.  Since  B(s)  is  in  the  denominator  of  each  response  factor 
expression,  we  must  find  the  roots  of  B(s)  =  0  in  order  to  find  the 
required  poles.  The  problem  is  complicated  even  for  a  two-layered 
si ab  for  which 

R?  , _  r _ _  R,  _ _  , _ 

B(s)  =  -  cosh  /sR7c7  sinh  /sR„C0  + - —  cosh  v'sR^Co  sinh  /sR, C, 

✓sR^  11  11  /  R1C1  2  2  11 

[Eq  2-76] 

We  are  quickly  forced  to  rely  on  numerical  techniques  to  find  the 
roots  of  B(s).  All  of  these  techniques,  however,  require  that  we  try 
various  values  of  s  in  the  appropriate  expression  for  B(s)  until  we 
find  s  such  that  B(s)  =  0.  Each  time  we  change  the  value  of  s  we 
must  perform  the  matrix  multiplication  necessary  to  calculate  each 
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element  of  the  transmission  matrix  (see  Equation  2-2 t).  Furthermore 
(theoretical ly) ,  an  infinite  number  of  roots  must  be  found.  (As  a 
practical  matter,  we  may  need  to  find  20  or  more,  depending  on  the 
properties  of  the  layers  in  the  slab.)  Even  after  the  roots  are 
found,  we  must  find  the  derivative  of  B(s)  with  respect  to  s 
evaluated  at  each  of  the  roots  (analogously  to  Equations  2-42  and 
2-54).  for  a  multilayered  slab,  we  must  carry  out  a  series  of  matrix 
man i pul ations  for  each  root  as  we  apply  t he  chain  rule  to  find  the 
required  derivatives.  Finally,  for  each  mth  response  factor,  we  must 
carry  out  the  sum  of  the  exponential  series  illustrated  by  the 
right-hand  sides  of  Equations  2-59  through  2-64.  Again,  in  theory, 
there  is  an  infinite  number  of  X,  Y,  and  Z  response  factors.  In 
practice,  20  or  more  response  factors  may  be  required  to  calculate 
heat  flow  through  heavy  masonry  walls  accurately. 

As  a  practical  matter,  the  calculation  of  response  factors  for 
multilayered  slabs  requires  the  use  of  a  computer. 

We  now  conclude  our  discussion  of  response  factors  by  presenting 
the  general  formula  for  calculating  each  term  for  a  multilayered 
si  ab. 


First,  let  -3  be  the  nth  root  of  B(s)  =  0. 
The  X  response  factor  series  is: 


II  O'"!  P 
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[Eq  2-82] 


The  2  response  factor  series  is: 
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[Eq  2-85] 


s=-Br 


where  A'(s)  =  — 


The  derivatives  A'(s),  B'(s),  and  D‘(s)  are  found  by  applying 
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the  chain  rule  as  follows: 
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[Eq  2-86] 


a;(s)  b'(s) 


C'(s)  D-<s) 


Finally,  we  note  two  properties  of  response  factors  which  are 
important  in  their  calculation  and  use.  First,  the  sum  of  each  of 
the  X,  Y,  or  Z  series  is  equal  to  the  U  value  of  the  composite  slab 
(U  =  1/[R1  +  R2  +  ...Rn]).  The  steady-state  limiting  case 

requires  these  sums  to  be  the  U  value  as  seen  by  Equations  2-70  and 
2-75.  In  Equation  2-75,  for  example,  if  both  T .  and  TQ  are  constant 
with  time,  then 


ii  t'-i  s 
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^0,1”  Yi  To  ^  Zm  [Eq  2-87] 

m=  1  m=  i 

Hence , 

l  Ym  =  f  Zm  =  U  LEq  2-88] 

m=  I  m= 1 

Similarly,  from  Equation  2-70, 

X  =  U.  [Eq  2-89] 

1  m 

This  has  proved  to  be  a  useful  check  in  calculating  response  factors. 

The  second  useful  property  of  response  factors  is  that  later 
terms  in  each  series  are  made  up  exclusively  of  exponential  functions 
(see  Equations  2-79,  2-82,  and  2-85).  Furthermore,  each  of  these 
terms  is  of  diminishing  importance,  since  the  roots  of  B(s)  =  0  are 
sequentially  more  negative.  We  will  make  immediate  use  of  this  pro¬ 
perty  in  defining  conduction  transfer  functions. 

Conduction  Transfer  Functions 

If  we  now  carefully  examine  the  response  factor  equations,  we 
see  that  the  higher-order  terms  have  the  same  basic  form.  For  the  Y 
response  factor,  for  example: 
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[Eq  2-90] 


for  m  >  3. 


on  the  negative  real  axis  and  which  increase  in  absolute  value  as  n 


increases.  For  a  heavy  single-layer  slab 
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Notice  that  for  m  =  20,  gj  A^0,  =  .0575,  but  g^0  =  -.00018  and 
g.A^^  -  U.OOOU.  Hence,  for  in  20,  Ym  -  g^A™.  Similarly,  for  m  = 
10,  Qj A J°  =  .2030,  g2A^°  =  -0.274,  but  g^0  =  .00023  and  g^J0  = 
0.0000.  Hence,  for  10  m  20,  Ym  is  approximately  equal  to  the 
contribution  caused  by  the  two  largest  roots,  i.e.,  Ym  -  ^  +  g^A^,1. 

The  general  pattern  is  clear.  For  large  m,  only  A^  is  important.  As 
m  decreases,  A^,  then  A^,  then  A^  and  so  on,  become  important.  A 
takes  its  value  from  the  root  of  B(s)  =  0  which  is  nearest  the  ori¬ 
gin,  A2  from  the  next  nearest  root,  and  so  on.  While  our  example  is 
for  a  single-layer  slab,  it  is  but  a  special  case  of  this  general 
behavior  of  later  terms  in  the  response  factor  series  for  multilay¬ 
ered  slabs. 


We  now  make  use  of  above  expressions  to  reduce  the  number  of 
response  factors  needed  to  calculate  heat  flux.  Consider  the  case 
where  the  heat  flow  at  the  inside  surface  is  to  be  calculated  under 
varying  outside  surface  temperature,  while  the  inside  surface  is  at 
zero  temperature. 


M  ,t 


m=l 


Y  T  .  xl 
m  o,t-m+l 


[Eq  2-92] 


If,  for  m  greater  than  or  equal  to  some  m',  Ym  -  g^A™,  we  may  write, 


for  times  t  and  t-1 : 
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[Eq  2-94] 


We  can  multiply  qi  t_1  by  and  subtract  the  result  from  q-  t. 
After  collecting  coefficients,  we  have: 

qi,t'Xlqi,t-l  =  *  YlTo,t  '  (Y2"X1  VTo,t-l  "  (Y3'XlY2)To,t-2  *••• 

'  (Ym,-rXlYtn,-2)To,t-m,+2  '  ^Ym,_XlYm'-l  )To,t-m'+l 

_(gl X!  "  91X1X1  ^To,t-m'  '  (glXl  _91X1X1  )To,t-m'-l 
-  ....  [Eq  2-95] 


Notice  that  the  coefficients  for  all  temperatures  preceding  T  ,  t-m+1 

{Ue"  To,t-m'  ’  To,t-m'-l  ’  To,t-m'-2’”  ^  are  zer0‘  Wc  may  therefore 
define  a  finite  series  of  first-order  conduction  transfer  coeffi¬ 


cients  or  first-order  conduction  transfer  function  as  follows: 
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Y 


1,1 


[Eq  2-96] 


I  ■} 

1  ,m 


=  Y  - 
m 


X,  Y 
1  m- 


for  in 


[Eq  2-97] 


We  may  now  rearrange  Equation  2-95  as: 
m' 


q. 


i,t  m^Yl  ,mTo, t-m+1  +  Xlqi,t 


[Eq  2-98] 


Applying  the  same  procedure  for  X  and  Z  response  factors  yields 
first-order  X  and  Z  conduction  transfer  functions: 


xl,l 

=  X1 

[Eq 

2-99] 

Xl,m 

=  Xm"XlVl 

for  m  >  2 

[Eq 

2-100] 

Zl.l 

=  Z1 

[Eq 

2-101] 

Z1  ,m 

=  Zm  '  XlZrn- 

1  for  m  >  2 

[Eq 

2-102] 

The 

counterparts  to  Equations  2 

-70 

and  2- 

75  then  become: 

m' 

m' 

°*i  ,t 

=  Y  X,  T. 
m=l  '’m’ 

, t-m+1  ’  J, 

m=  I 

Y 

T 

m  o 

,  t-m+1 

[Eq 

2-103] 

m1 

m' 

qo,t 

=  y  Y,  T. 

Li  1  ,m  l 
m=l 

"  m=l 

Y 

T 

m  o 

,t-m+l 

[Eq 

2-104] 

Aj  is  often  termed  the  "common  ratio"  since,  for  large  m,  it 
equals  the  ratio  of  successive  terms  in  the  response  factor  series 
(i.e.,  Y^j/Y^).  It  is  neither  necessary  nor  advisable,  however,  to 


calculate  A^  as  the  ratio  of  successive  terms,  since  A^  is  known  from 
calculations  which  precede  the  final  step  used  to  calculate  response 
factors.  It  is  also  not  necessary  to  decide,  a  priori,  on  a  value 
for  m'.  Instead,  a  reasonable  procedure  is  to  calculate  response 
factors  until  their  values  are  some  very  small  function  of  the 
overall  U  value  for  the  multilayered  slab.  Conduction  transfer  func¬ 
tions  can  then  be  calculated  and  the  series  may  be  truncated  when  the 
mth  conduction  transfer  function  is  some  very  small  fraction  of  the  U 
value  times  (1  -A^).  Figure  2-7  shows  a  plot  of  the  results  of  such 
a  procedure  for  the  Y  conduction  transfer  function,  using  the  heavy, 
single-layer  slab  of  our  previous  example.  The  upper  solid  line 
shows  the  Y  response  factors,  and  the  upper  dashed  line  shows  the 
first-order  conduction  transfer  function. 

The  lower  solid  and  dashed  lines  of  Figure  2-7  are  second-  and 
third-order  conduction  transfer  functions,  a  logical  extension  of 
first-order  conduction  transfer  functions.  In  going  from  response 
factors  to  first-order  conduction  transfer  functions,  we  have 
accounted  for  the  effects  of  all  response  factors  from  m1  to  infinity 
in  one  flux  history  term  and  removed  the  effects  of  A^  from  the 
remaining  conduction  transfer  functions.  We  are  left  with  first- 
order  conduction  transfer  functions  whose  later  terms  (as  m  -*  m') 
are  now  dominated  by  A,,.  We  now  remove  the  effect  of  A^>  in  a  com¬ 
pletely  analogous  way.  The  necessary  sample  equation,  similar  to 
Equation  2-95,  is: 


Figure  2-7.  Y  concJuct.ion  transfer  functions  for  single-layer  slab 
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(qi,t'Alqi,t-l^  '  X2^qi ,t-l _Xlqi ,t-2)  ~ 


-Y,  ,T  -  (Y,  9-X„Y,  i )T  ,  ,  -  (Y.  _-A0Y,  ,)T  .  , 

1,1  0 , t  1,2  2  1,1  o , t- 1  1,3  2  1,2  o,t-l 


‘  (Yl,m"'A2Yl,m',-l)To,t-in,,+l 


[Eq  2-105] 


where  m' 1  is  the  point  after  which  Y,  „  ~  gJ?.  Second-order  conduc- 
tion  transfer  functions  can  now  be  defined  completely  analogously  to 
first-order,  conduction  transfer  functions,  as  can  third-,  fourth-, 
or  kth-order  conduction  transfer  functions. 


X2,l 

=  xi 

,i 

[-Q 

2-106] 

X2,m 

=  xi 

-A0X-.  m  i  for  m  >  2 
,m  2  1  ,m- 1  — 

[Eq 

2-107] 

Y2 ,1 

=  Yi 

,1 

[Eq 

2-108] 

Y2,m 

=  Yi 

m  ■  A0Y1  mi  for  m  >  2 
,m  2  1  ,m- 1  — 

[Eq 

2-109] 

r>i 

no 

=  zi 

,1 

[Eq 

2-110] 

Z2,m  Z1 ,m  "  X2Z1 ,m-l 


for  m  >  2 


[Eq  2-111] 
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The  counterparts  to  Equations  2-103  and  2-104  for  second-order  con¬ 
duction  transfer  functions  are: 

m"  m" 


^2  ,nJi  ,t-m+l 
m=  1 

m" 

^  Y2,niT2,t-m+l 
m=  I 


T  Y„  T  .  Al 
^1  2 ,m  o,t-m+l 

m" 

n!=1  ^2,m^o,t-m+l 


+  ^xi+X2^qo,t-TXlqo,t-2 


[Eq  2 


[Eq 


For  kth-order  conduction  transfer  functions,  we  have 


ri  ii 

q1,t  "  m|i Xk ,mTi  ,t-n+l  ’  J;-,  Yk,mTo,t-m+l  +  Jj/n^i.t-k 


m=l 


m=l 


[Eq  i 


qo,t  Yk,inTi  ,t-m+l  ^ Zk,mTo,t-m+l  +  m^Fmqi,t-k 


[Eq 


where  M  is  the  value  of  m,  above  which  the  kth  conduction  transfer 
function  is  proportional  to  A,171.  F  values  are  flux  history  coeffi- 

l\  ill 

cients  defined  as  follows: 
k 

F,  =  I  Xn 
1  n=l  n 


-112] 


-113] 


-114] 

2-115] 


[Eq  2-116] 


F2  =  -  [sum  of  the  products  of  A taken  two  at  a  time] 


[Eq  2-117] 


=  +  [sum  of  the  products  of  A^  taken  three  at  a  time]  [Eq  2-118] 

k+l  ^ 

Fk  =  (-1 r  1  n  A  [Eq  2-1131 

K  n=l  n 

For  k  =  4,  for  example  (fourth-order  conduction  transfer  functions): 

F1  =  A1  +  A2  +  A3  +  A4  [Eq  2-120] 

f2  ~  '  +  A1A3  +  A1A4  +  A2A3  +  A2X4  +  A3A4^  [F<1 

F3  =  +  [A-,A2A3  +  A1A2A4  +  A1A3A4  +  A2A3A4]  [Eq  2-122] 

F4  =  •  A1A2A3A4  [Eq  2_123] 


The  practical  value  of  conduction  transfer  function  and  a  gen¬ 
eral  scheme  for  determining  the  order  of  the  transfer  function  to  use 


are  both  illustrated  by  Figure  2-7.  First,  we  see  that  by  using  one 
flux  history  term,  we  can  use  about  20  temperature  history  terms 
instead  of  about  50.  Thus,  to  calculate  fluxes  using  first-order 
conduction  transfer  functions,  we  must  perform  about  86  multiplica¬ 
tions  and  additions  and  keep  track  of  about  40  temperature  and  flux 
histories,  compared  to  about  200  multiplications  and  additions  and 
about  98  temperature  histories  if  we  had  used  response  factors.  For 
second-order  transfer  functions,  about  10  past  temperatures  are 
required,  along  with  two  past  fluxes  (46  operations  and  22  his¬ 
tories).  For  third-order,  about  eight  past  temperatures  are  needed, 
along  with  three  past  fluxes  (42  operation  and  20  histories).  Obvi¬ 
ously,  the  reduction  in  the  number  of  terms  required  is  diminishing 
at  this  point,  but  we  have  already  reduced  the  needed  calculations 
and  storage  to  about  20  percent  of  what  would  be  required  if  we  had 
used  response  factors  directly.  Our  efforts  would  be  half  of  what 
would  be  required  with  first-order  conduction  transfer  functions. 

If,  in  the  example  shown  in  Figure  2-7,  we  proceed  to  fourth-  or 
fifth-  or  higher-order  transfer  functions,  we  will  first  simply  trade 
one  pair  of  past  temperatures  for  a  pair  of  past  fluxes.  Soon,  how¬ 
ever,  the  number  of  flux  terms  required  will  increase  more  rapidly 
than  the  reduction  in  the  number  of  required  temperature  histories. 

In  the  limit,  we  will  need  an  infinite  number  of  flux  terms,  along 
with  the  current  temperature.  The  selection  of  the  appropriate 
number  and  order  of  conduction  transfer  function  coefficients  to  use 

j  j 

i 


II  II 


for  a  given  wall  section  has  been  addressed  by  Mitalis  (ref  24),  and 
by  Peavy  (ref  18)  and  Walton  (ref  25).  (The  mathematical  development 
of  the  procedures  used  by  Mitalis,  Kusuda  [ref  17],  and  Peavy  are 
fundamentally  identical.) 


To  determine  the  appropriate  number  of  terms  for  kth-order  con¬ 
duction  transfer  functions,  we  recall  that  for  zero-order  conduction 
transfer  functions  (i.e.,  response  factors),  the  sum  of  the  coeffi¬ 
cients  of  the  infinite  X,  Y,  or  Z  series  equals  the  U  value.  For 
first-order  conduction  transfer  functions. 


1  ,m 


Mi 

“  I  Y 

m=l  1 


-  u 


m=l 


1  ,m 


(1-A-,)U 


[Eq  2-124] 


(approximately).  Similarly,  for  second-order  conduction  transfer 
functions , 


2  ,m 


M2 

=  l  Y 


i  2  ,m 

m=l 


'  I  z2.m  =  <>-x1)<l-*2)u 


m=l 


[Eq  2-125] 


In  general , 


[Eq  2- 
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l  Y 
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k  ,tn 
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I  Y 
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=  l 
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'k  ,'ii 


=  U 


k 

[I 

i=1 


(1-Ai) 


M  can  be  determined  by  requiring  that  the  sum  of  the  first  M  terms  of 
each  of  the  X,  Y,  or  Z  conduction  transfer  functions  be  equal  to  a 
large  fraction  (say  99  percent)  of  the  right  side  of  Equation  2-125. 

Two  approaches  for  determining  k  --  the  order  of  the  transfer 
functions  --  have  been  proposed.  Mitalis  suggests  that  the  order  of 
the  conduction  transfer  functions  should  be  one  less  than  M  (i.e., 
k  =  M-1).*  Extensive  tabulations  of  conduction  transfer  functions  of 
order  k  =  M-1  can  be  found  in  the  ASHRAE  Handbook  of  Fundamentals 
(ref  21).  Peavy  (ref  18)  suggests  that  k  need  be  no  larger  than  5 
for  common  (even  massive)  walls,  roofs,  and  floors.  Using  5  as  an 
upper  limit,  Walton  (ref  25)  suggest.,  that  first  values  of  X,  Y,  and 
Z  conduction  transfer  function  coefficients  be  calculated  for  all  six 
possible  orders  (i.e.,  M  =  1,  k  =  0,  1,  2,  3,  4,  and  5;  then  M  =  2, 
k  =  0,  1,  2,  3,  4,  and  5,  and  so  on).  Each  time  M  is  incremented, 
the  test  of  Equation  2-125  is  applied  for  successively  increasing 

*  This  approach  is  a  natural  consequence  of  Mitalis1  use  of  the  Z 
transform  (not  to  be  confused  with  the  Z  conduction  transfer  func¬ 
tion)  and  polynomial  arithmetic.  It  does  not  necessarily  lead  to  a 
minimum  number  of  coefficients  for  a  given  precision  requirement. 
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values  of  k.  (As  a  refinement,  we  note  that  the  test  need  only  be 
made  for  k  <  M.)  As  soon  as  Equation  2-125  is  satisfied  to  the 
required  precision,  k  and  M  are  both  known.  Figure  2-8  shows  a  flow 
diagram  of  the  process.  For  a  given  precision,  Walton's  procedure 
yields  the  smallest  possible  set  of  conduction  transfer  functions 
and,  therefore,  results  in  the  minimum  number  of  calculations  when 
flux  calculations  are  made  using  the  conduction  transfer  functions. 


CALCULATE 
Xk,M*  Yk,M’Zk,M 


3  FREQUENCY  RESPONSE  OF  MULTILAYERED  SLABS 


Exact  Solution  From  the  Transmission  Matrix 

We  now  examine  the  special  case  of  a  wall  with  fixed  (zero) 
inside  surface  temperature  and  sinusoidally  varying  outside  tem¬ 
peratures.  This  case  is  examined  for  two  reasons.  First,  a  large 
portion  of  the  variation  in  ambient  outdoor  temperature  and  solar 
radiation  is  periodic  and  sinusoidal  in  character.  These  two  cli¬ 
mate  variables  can  be  used,  along  with  the  absorptivity  of  the  out¬ 
side  surface  of  a  wall  or  roof,  to  establish  the  so-called  sol -air 
temperature.  If  we  treat  the  thermal  resistance  of  the  air  film  ad¬ 
joining  the  outside  surface  of  a  multilayered  slab  as  another  layer 
of  the  slab,  then  the  sol -air  temperature  becomes  the  "outside  sur¬ 
face"  boundary  condition.  The  annual  variation  in  the  sol -air  tempera¬ 
ture  for  a  given  location  and  exposure  can  be  approximated  as  the 
sum  of  a  relatively  few  sine  and  cosine  terms  at  various  fixed  fre¬ 
quencies.  This  simplicity  is  in  sharp  contrast  to  the  use  of  hourly 
discrete  time  series  containing  8760  hourly  values  per  year. 

The  second  reason  for  examining  this  case  is  that  the  special 
instance  of  one  sinusoidal  and  one  fixed  temperature  boundary  condition 
allows  us  to  solve  the  heat  conduction  equation  much  more  easily,  as 
will  be  demonstrated  in  subsequent  paragraphs. 

Recall  the  general  matrix  form  of  the  Laplace  transform  solution 
of  the  neat  conduction  equation  for  multilayered  slabs  [Equation  2-27]: 
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T-j  (s)~ 

A(s) 

B{  s  )” 

'Ws>' 

_^1 (s) 

C(s) 

D(s)_ 

-Vl(s)- 

[Eq  3-1] 


Since  we  are  interested  in  the  resulting  fluxes  when  temperature 
boundary  conditions  are  specified,  we  rewrite  Equation  3-1  as  follows: 


q-,  (s) 

-Vi(s>- 


r  dui 

B(sT 

i 

L  bTsT 


bTsT 


~T,(s) 

-Ws>- 

[Eq  3-2] 


For  the  special  case  of  T-|(s)  =  0,  we  may  find  q-j(s)  as: 


[Eq  3-3] 


The  transfer  function  for  the  wall  becomes: 
q,(s)  , 

Vi111  '  '  ^ 


[Eq  3-4] 


We  now  wish  to  find  q ^ ( s )  under  steady  periodic  state  when 

Tn+-j  ( t )  =  sin(wt)  [Eq  3-5] 

where  w  =  angular  frequency  of  temperature  variation  of  surface  T  +i . 

A  well-known  property  of  the  Laplace  transform  is  that  if  s  =  ju>,  the 
Laplace  transform  becomes  the  Fourier  transform.  Under  the  special 
case  that  Tn+j(t)  =  sin(wt)(for  -  «  t  +  °°)  and  T-j(t)  =  0.0, 
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the  inverse  transform  of  Equation  3-3  is: 
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the  magnitude  of  the  complex  valued  function  - 


=  the  argument  of  -  g^j^y  • 


[Eq  3-6] 


[Eq  3-7] 
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Thus,  q-j(t)  can  be  viewed  as  a  weighted  sum  of  sin(cot)  and  cos(wt)  or 
as  a  sine  wave  of  amplitude  |-  gf j^y  I  and  phase  lag  <j>.  We  can  write 
q-|(t)  for  a  unit  sinusoidal  outside  surface  temperature  variation  of 
frequency  w  as  either: 
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[Eq  3-8] 


q-|(t)  =  U^sin(wt-<})) 


[Eq  3-9] 


where  U  =  /W  +  U 
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We  immediately  see  how  much  easier  it  is  to  calculate  heat  flux 
in  this  special  case.  To  find  U  „  and  U  .  we  simply  substitute  jw 
for  s  in  each  term  of  the  transmission  matrix  expressions  for  each  • 
layer,  carrying  out  the  complex  matrix  multiplications  to  find  B(jw), 
and  find  the  real  and  imaginary  parts  of  -  'g^y  •  for  example,  for  a 
two-layered  slab: 


A(j<>0  B(  j(»>j 

_C(joj)  D(jw] 


_ Rt  _ 

cosh/jwR,  C,  ,  - -  sinh/j<i)R,C, 
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\  1  V\/ju,R2C2  c 

+  ( — 1  —  si  nhyjujRTcTj  (  cosh/jwR9C, 
V/jcuR1C1  1  '/  V  1  1 


[Eq  3 


Expressions  for  slabs  of  more  than  two  layers  become  complicated,  but 
B(jw)  can  easily  be  calculated  with  simple  computer  algorithms. 
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Compared  to  the  response  factor  calculation,  the  frequency  response 
of  multilayered  slabs  can  be  calculated  without  (1)  numerically 
searching  for  the  roots  of  B(s)  =  0,  (2)  finding  derivatives  of  B(s) 
at  the  roots  of  B(s)  =  0,  and  (3)  evaluating  and  summing  residues  at 
B(s )  =  0. 

Figure  3-1  shows  U  ,  the  frequency-dependent  conductance  or  U- 

value,  and  the  phase  lag,  <j>,  versus  frequency  for  a  heavy,  single- 

m^°K 

layer  slab  with  inside  and  outside  air  film  resistances  of  .12  - - 

o  W 

rn  °K 

and  .029  — -,  respectively  (essentially  three-layered  slabs). 

Figure  3-2  shows  similar  plots  for  a  heavy,  multilayered  wall  section 
with  insulation  separating  the  heavy  layers. 

Equivalent  Single-Layer  Thermal  Properties 

One  of  the  objectives  of  this  research  was  to  determine  whether 
multilayered  slabs  could  be  approximated  as  single-layer  homogeneous 
slabs  having  "equivalent"  thermal  properties.  To  this  end,  the  error 
associated  with  such  an  approximation  was  studied.  There  were  two 
reasons  for  considering  multilayered  slabs  as  having  equivalent  single¬ 
layer  resistance,  Rg,  and  equivalent  single-layer  capacitance,  Cg: 

(1)  the  potential  for  simplifying  the  calculation  of  B(jw),  and  (2) 
the  potential  for  experimentally  determining  values  for  C£  and  Rg  for 
complicated  wall  sections. 

We  now  consider  a  two-layered  slab  and  define  Rg  and  Cg  for  this 
case,  arguing  that  if  these  equivalent  parameters  can  be  found  for  a 
two-layered  slab,  the  procedure  can  be  recursively  applied  to  find  Rg 
and  Cg  for  n  layers. 


PHASE  LAG  MAGNITUDE 


FREQUENCY,  CYCLES  PER  DAY 


Returning  to  the  basic  definition  of  the  terms  in  the  trans¬ 


mission  matrix,  we  have  for  a  single-layer  slab: 


B(s)  =  ■  sinh  /sRTC 

/sR-j  C-j 


For  a  two-layered  slab,  we  have: 


B(s )  =  A-|  ( s ) B2 ( s )  +  B1(s)D2(s) 


_  R?  _ 

=  (cosh  ZsRtC,  )!  - sinh  /sR0C„ 

V  1  V  VsR2c2  j 

+  ( - —  sinh  /sR,C,  ]  (cosh  AR0C, 

\/sR, C,  V  V 


Now,  if  equivalent  single-layer  thermal  properties,  Rg  and 
Cg,  exist  for  a  two-layered  slab,  then,  from  Equations  3-12  and 
3-13,  they  clearly  must  satisfy  the  following  condition: 


Unfortunately,  while  we  can  solve  Equation  3-14  for  Rg  and  Cg  by 
equating  real  and  imaginary  parts  to  form  two  equations  in  two  un¬ 
knowns,  we  note  that  R  and  C  are  not  independent  of  s.  (This 

e  e 


[Eq  3- 


[Eq  3- 


[Eq  3 


79 


implies,  since  s  =  jw  in  our  previous  expression  for  q(t)  for  the 
sinusoidal  exterior  boundary  condition,  that  Rg  and  are  frequency- 
dependent.)  The  Taylor  series  expansion  of  Equation  3-14  will 
clarify  this: 


Let  Zi  =  /sFL  C.  - 


[Eq  3-15] 


The  series  expansion  of  sinh  Z  and  cosh  Z  are  as  follows: 


Z^  Z^  I? 
sinh  z  =  Z  +  j r  +  5T  +  7T 


Z^  z4  Z^ 

cosh  Z  =  1  +  21  +  4j  +  61 


[Eq  3- 16] 

[Eq  3-17] 


Using  Equations  3-15,  3-16  and  3-17  in  Equation  3-14  provides: 


7  +  _s.  + 

e  3!  5!  •' 
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z3  ^5 

z  +  -1  +  -I 

L]  3!  5! 


Z2  Z4  ' 

1  +  11+  II 
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[Eq  3-18] 


Rearranging  and  simplifying  produces: 


2  3  2  2 

Re  s  RJcV+... 

r  +  e.  e-  + 

Ke  3!  5! 


=  ( R-j  +Rg )  + 


(r2C1  +3R1 R2C1  +3R-|  R2C2+R2C2)s 
_ 


(R2C2+5R1 R2C2+5R2R2C2+R2C2)s2+. 
— 


[Eq  3-19] 


We  now  see  that  if  we  approximate  Equation  3-19  by  using  only  the 
first  two  terms  of  both  sides,  we  can  define  Rg  and  Cg  independently 
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of  s.  Thus,  we  can  define  R  and  C  as  constant  equivalent  thermo- 

e  e 

physical  properties  only  in  this  approximate  sense. 

We  now  consider  whether,  and  under  what  circumstances.  Equation 

3-19  can  be  approximated  by  using  only  the  first  two  terms.  Under 

this  approximation,  R  and  C  are  defined  as: 

e  6 

Re  =  R-|  +  r2  [Eq  3-20] 

R?C,  +  3R, R9C,  +  3R,R9C9  +  R^C„ 

C  =  -!-! - — — o— -L-  - - -  [Eq  3-21] 

Re 

We  can  apply  Equations  3-20  and  3-21  recursively  for  slabs  with  more 
than  two  layers  by  first  calculating  Rg  and  Cp  for  the  two  outermost 

layers;  then,  treating  the  outer  two  layers  as  one  equivalent  layer, 

we  substitute  Rg  for  R-j ,  Cg  for  C-| ,  and  R3  and  C3  (properties  of  the 
third  outermost  layer)  for  R3  and  and  calculate  new  values  for  Rg 

and  Cg,  etc.  Since  we  can  calculate  Rg  and  Cg  for  any  multilayered 

slab,  we  can  compare  the  frequency  response  of  a  multilayered  slab, 
calculated  by  the  equivalent  single-layer  slab  model,  to  the  exact 
solution,  obtained  by  evaluating  the  terms  of  the  transmission  matrix. 

To  compare  them,  we  set  s  equal  to  jw,  and  plot  |-  anc*  $  for 

various  frequencies,  using  both  models.  These  comparisons  are  shown 
for  several  multilayered  slabs  of  practical  interest  in  Figures  3-3 
through  3-9. 

Figure  3-3  shows  results  for  a  single  10.16-cm  (4-in.)  con¬ 
crete  slab  with  outside  and  inside  air  film/surface  resistances  of 

2  2 
m  °  K  m  °  k 

.029  — ^  and  .12  ,  respectively.  (These  air  film  resistances 
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Figure  3-3.  Frequency  response  of  heavy  concrete  slab. 
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Figure  3-6.  Frequency  response  of  a  multilayered  slab, 
brick  and  heavy  concrete. 
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are  combined  radiative-convective  resistances  and  are  used  for  all 
wall  sections  studied.)  We  see  that  the  single-layer  equivalent  slab 
model  agrees  reasonably  well  with  the  multilayered  model.  As  ex¬ 
pected,  the  agreement  at  low  frequencies  is  nearly  exact,  since  it  is 
in  this  range  (small  jw)  that  the  two-term  approximations  for  sinh  and 
cosh  are  very  close  to  the  actual  values.  However,  we  have  little 
interest  in  this  region,  since  it  represents  what  is  essentially 
steady-state  heat  transfer.  As  shown  in  Chapter  4,  variation  in 
climate  is  characterized  by  periodic  variations  at  daily  and  twice- 
daily  frequencies.  Therefore,  if  we  are  to  use  the  single-layer  equi¬ 
valent  slab  model,  its  agreement  with  the  multilayered  slab  model  must 
be  usable  at  frequencies  of  one  cycle  per  day  and  two  cycles  per  day. 
For  the  slab  in  Figure  3-3,  the  single-layer  model  magnitude  is  in 
error  by  about  6  percent  at  one  cycle  per  day,  and  about  18  percent  at 
two  cycles  per  day. 

Much  worse  agreement  is  shown  in  Figure  3-4,  which  compares  the 
single-layer  slab  model  to  the  exact  solution;  this  comparison  is  for  a 
multilayered  slab  consisting  of  10.6-cm  face  brick,  7.62-cm  insulation, 
and  10.16-cm  heavy  concrete.  In  this  case,  the  single-layer  slab  model 
predicts  a  magnitude  which  is  nearly  twice  that  of  the  exact  solution 
at  a  frequency  of  one  and  two  cycles  per  day.  The  phase  lag  is  also 
in  error  by  the  equivalent  of  about  3  hours  at  one  cycle  per  day 
and  1.5  hours  at  two  cycles  per  day. 

It  was  found  that  consistent,  and  often  considerable,  improvement 
in  agreement  between  the  equivalent  single-layer  model  and  the  exact 


model  could  be  achieved  by  combining  the  inside  air  film/surface  re¬ 
sistance  with  the  resistance  of  the  Innermost  wall  layer.  Calculating 
the  equivalent  capacitance  in  this  way  for  the  wall  section  of 
Figure  3-4  yields  the  much  improved  results  shown  in  Figure  3-5. 

Figures  3-6  through  3-8  compare  several  common  wall  sections  for 
the  modified  equivalent  single-layer  slab  model  and  the  exact  model. 
Figure  3-6  illustrates  a  wall  section  made  of  two  different,  but 
relatively  heavy,  masonry  materials.  Figure  3-7  shows  a  heavy  wall 
layer  that  is  exposed  on  the  outside  but  insulated  on  the  inside. 

Figure  3-8  shows  a  lightweight  wall  section  (low  on  "thermal"  mass), 
and  confirms  that  heat  flow  through  light  walls  can  be  analyzed  using 
steady-state  procedures;  the  magnitude  of  B(jco)  is  nearly  constant,  and 
the  phase  lag  is  nearly  zero. 

Figure  3-9  shows  that  even  with  the  modification  used  to  calculate 
Cg,  severe  errors  can  occur  when  the  equivalent  single-layer  model  is 
used.  This  figure  shows  typical  results  obtained  when  a  heavy,  con¬ 
ductive  layer  is  insulated  on  the  outside.  This  type  of  wall  construc¬ 
tion,  though  not  common  in  existing  buildings,  may  prove  important  in 
future  energy-conservative  designs. 

Figures  3-3  through  3-9  are  but  samples  of  the  many  comparisons 
made  between  the  equivalent  single-layer  model  and  the  exact  model  when 
calculating  the  frequency  response  of  multilayered  slabs.  However, 
they  provide  enough  information  to  evaluate  the  equivalent  single-layer 
model . 

It  is  therefore  concluded  that  the  equivalent  single-layer  model 
is  a  potentially  dangerous  means  of  quantitatively  estimating  the 
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dynamic  behavior  of  multilayered  slabs.  For  important  frequencies 
and  for  realistic  wall  sections,  the  arguments  of  the  hyperbolic 
functions  making  up  the  transmission  matrix  are  simply  too  large  to 
be  approximated  by  a  two-term  expansion.  For  the  same  reason,  the 
frequency  response  of  multilayered  slabs  cannot  be  characterized  with 
reasonable  accuracy  by  any  model  simpler  than  the  exact  model  used 
as  the  standard  for  comparison.  Furthermore,  the  inclusion  of  more 
terms  in  the  series  expansions  leading  to  the  equivalent  single-layer 
mode1  forces  us  to  acknowledge  the  true  frequency  dependence  of  the 
equivalent  single-layer  thermodynamic  properties.  Calculating  fre¬ 
quency-dependent  values  for  Rg  and  Cg  involves  the  solution  of  simul¬ 
taneous  non-linear  equations  (like  Equation  3-14)  and  is  more  laborious 
than  calculating  the  exact  solution  itself. 

The  first  section  of  this  chapter  showed  how  the  frequency  re¬ 
sponse  of  a  multilayered  slab  can  be  calculated  much  more  easily  than 
response  factors.  In  this  second  section  we  have  shown  that  attempts 
to  further  simplify  the  dynamic  heat  conduction  problem  were  only  of 
limited  value,  and  furthermore,  given  the  simplicity  with  which 
frequency  response  can  be  calculated  from  the  transmission  matrix, 
further  simplifications  are  probably  unwarranted. 
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4  PERIODIC  AND  STOCHASTIC  BEHAVIOR  OF  WEATHER  DATA 

Introduction 

In  Chapter  3,  the  calculation  of  heat  conduction  through  multi¬ 
layered  building  walls  and  roofs  was  shown  to  be  much  less  compli¬ 
cated  if  the  outside  surface  temperature  (or  the  sol -air  temperature) 
is  treated  as  the  sum  of  pure  sinusoids.  This  requires  that  the 
actual  outside  temperature  and  solar  radiation  be  adequately  approxi¬ 
mated  by  a  small  set  of  sine  and  cosine  functions  oi  time.  In  this 
chapter,  we  analyze  actual  weather  data  to  determine  whether  or  not 
such  an  approximation  is  reasonable.  Ihe  following  sections 
describe:  (1)  the  weather  data  chosen  for  analysis,  (2)  three 
analysis  techniques  applied  to  characterize  weather  behavior,  and  (3) 
the  results  of  using  these  techniques. 

Weather  Data 

Hourly  weather  data  for  four  different  locations  in  the  United 
States  were  obtained  from  the  Environmental  Technical  Applications 
Center  (ETAC)  Air  Weather  Service,  United  States  Air  Force.  The 
tapes  selected  were  of  the  so-called  SOLMLT  format.  These  tapes  con¬ 
tain  about  20  years  of  measured  hourly  climate  data.  The  four  sites 
selected  represented  four  distinctly  different  climates:  Charleston, 
SC,  represented  a  hot  humid  climate;  Madison,  WS,  a  cold  northern 


climate;  Fort  Worth,  TX,  a  hot  dry  climate;  and  Santa  Maria,  CA,  a 
marine  cl imate. 


From  the  20  years  of  data  available  for  each  site,  a  lC-year 
period  covering  the  years  1 9E>3  through  196k?  for  Charleston,  Madison, 
and  Fort  Worth,  and  ttie  years  1954  through  1963  for  Santa  Maria  was 
chosen  tor  analysis.  Three  types  of  data  time  series  were  analyzed 
to  investigate  short-  and  long-term  data  trends. 

1.  Thirty-day  periods  of  hourly  data  and  40-day  periods  of  dat 
taken  every  other  hour 

2.  Ten  years  of  daily  average  data  (3650  samples  per  weather 
data  variable) 

3.  Ten  years  of  5-day  average  data  (730  samples  per  weather 
data  v  a  r i a  b 1 e ) 

The  two  climate  variables  that  affect  building  heat  transfer 
most  directly  (contribute  to  the  sol-air  temperature)  are  dry-bulb 
temperature  and  beam  solar  radiation.  Hence,  both  of  these  time 
series  were  analyzed. 

Procedures 

in  order  to  characterize  climate  data,  we  must  determine  whether 


the  data  is  the  result  of  a  determini  Stic  process  or  a  stochastic 


r 


process.  If  we  can  develop  a  model  which  enables  us  to  calculate  the 
value  of  a  time-dependent  variable  (dry-bulb  temperature,  for  exam¬ 
ple)  nearly  exactly  at  any  instant  of  time,  the  model  and  the  process 
ore  cal  led  determim  Stic:.  Box  and  Jenkins  (ref  2b)  use  the  example 
of  a  projectile  trajectory  which  can  be  calculated  nearly  exactly  if 
the  direction  and  velocity  at  launch  are  known.  A  stochastic  process 
is  one  which  cannot  be  modeled  deterministically  because  there  are 
unknown  factors  affecting  the  variable  which  prevent  the  exact  calcu¬ 
lation  of  its  future  behavior.  We  can  develop  stochastic  models  for 
these  processes  which  allow  us  to  calculate  the  expected  future  value 
of  the  process  and  the  probabi 1 i ty  that  a  future  value  will  be 
between  two  specified  limits. 

There  are  also  certain  time  series  which  can  be  viewed  as  being 
the  result  of  both  deterministic  and  stochastic  processes.  A  radio 
signal  containing  "static"  is  a  familiar  example.  We  expect  that 
weather  data  can  also  be  modeled  as  being  the  result  of  both  deter¬ 
ministic  and  stochastic  processes.  Extraterrestrial  (or  "clear  sky") 
radiation,  for  example,  can  be  calculated  from  basic  laws  of  astron¬ 
omy  (deterministically).  The  occasional  passing  of  clouds  cannot  be 
exactly  predicted,  however;  this  must  be  viewed  as  a  stochastic  pro¬ 


cess. 
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In  addition  to  our  a  priori  expectation  that  both  known  and  un¬ 
known  factors  contribute  to  climate  variation,  we  must  also  recognize 
that  we  are  searching  for  climate  models  which  are  to  be  used  for  the 
specific  purpose  of  predicting  energy  use.  If  we  are  to  take  advan¬ 
tage  of  the  results  described  in  Chapter  3,  we  need  a  model  which 
uses  sines  and  cosines  to  describe  climate  variation.  The  three  time 
series  analysis  techniques  described  below  were  chosen  with  this 
requirement  in  mind. 


The  first  approach  used  to  analyze  the  weather  data  time  series 
was  the  discrete  Fourier  transform,  which  is  one  of  the  most  direct 
methods  of  analyzing  discrete,  periodic  data.*  For  a  real  data 
sequence,  K-  of  N  equally  spaced  samples,  the  discrete  Fourier 
transform  of  X-  consists  of  N  coefficients  ar  and  br>  r=0,  1.2...N/2 
such,  that 


N/2-1 

X.  =  a  +  2  Y  a  cos(2TTri/N) 

1  0  r=  1  r 

N/2-1 

+  2  Y,  b  sin(2irri/N) 
r=  1  r 

+  a^  cos(2mj N/2N) 

*  The  direct  use  of  the  Fourier  transform  would  not  be  appropriate 
for  a  time  series  which  exhibited  only  stochastic  behavior  (ref 
27).  We  know,  however,  that  the  earth  revolves  about  its  axis  and 
orbits  the  sun  in  a  completely  predictable  way.  We  assume  that 
this  predictable  behavior  causes  part  of  the  variation  in  climate 
(the  deterministic  part).  We  can,  therefore,  expect  to  gain  mean¬ 
ingful  insight  into  the  behavior  of  the  weather  data  time  series, 
particularly  their  deterministic  components,  by  the  direct  use  of 
the  discrete  Fourier  transform. 


where  i  denotes  the  ith  sample  beginning  at  0.  For  a  series  sampled 
over  t i me ,  denotes  the  value  of  the  series  at  time  i’t,  where  t 
is  the  sampling  interval.  The  coefficients  ar  and  br  are  determined 
from  the  data  as  follows: 

,  N-l 

a  =  tt  z  X,  cos(2irri/N)  [Eq  4-2] 

1=0 

1 

br  =  ¥  Y  Xi  cos(27iri/N)  [Eq  4-3] 

r  "  i=0  1 

We  see  from  Equation  4-1  that  if  many  of  the  coefficients  a^  and 
br  turn  out  to  be  small  for  a  giver:  data  set,  then  X.  can  be  approxi¬ 
mated  by  taking  the  sum  of  a  few  sine  and  cosine  terms. 

To  compare  the  importance  of  the  various  terms  of  the  data  sets 
under  study,  the  discrete  Fourier  transform  was  calculated  by  using 
computer  subroutines  available  in  the  International  Mathematics  and 
Statistics  Library  (IMSL),  a  subprogram  library  available  at  most 
computer  centers.  Computer  software  was  developed  to  plot  the  magni¬ 
tude  of  the  contribution  of  a  and  b  at  each  frequency.  Specifically, 
aQ  (the  mean  of  the  data  set)  was  plotted  at  zero  frequency.  The 

quantity  2  >^r  +  b^  was  plotted  for  discrete  frequencies  between  one 
cycle  per  sampling  period  and  (N/2-1)  cycles  per  sampling  period. 
aN/2  was  plotted  for  the  Nyquist  frequency  of  N/2. 


Note  that  the  so-called  power  spectrum  in  its  discrete  form 
would  be  plotted  as  a^  +  b^.  The  decision  not  to  use  the  power  spec¬ 
trum  form  stems  from  the  physical  interpretation  which  can  be  given 
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to  the  quantities  aQ  and  2  /ar  +  b^.  aQ  is  the  mean  value  of  the 
variable  (°C  for  temperature,  w/m^  for  radiation).  The  quantity 

r?  t 

2  /a  +  b  is  the  amplitude  of  the  contribution  to  the  variation 
r  r 

about  the  mean  due  to  a  sinusoid  of  frequency  r  cycles  per  sampling 
period. 

The  second  analysis  approach  used  consisted  of  finding  adequate 
autoregressive  moving  average  (ARMA)  models  for  the  temperature  and 
beam  solar  radiation  time  series.  An  ARMA  model  of  order  m,n  (an 
ARMA  (m,n)  model)  has  the  form 


'lXt-l  + 


+  *  mXt-m 


+  at  '  Vt-l  "  92at-2 


0nat-n 


[Eq  4-4] 


where  is  the  value  of  the  discrete  time  series  variable  at  time  t, 
and  Xt  .  is  the  value  of  the  time  scries  at  time  t-  At,  where  At  is 
the  sampling  interval.  Xt  ^  is  the  value  of  the  time  series  variable 
at  t-2  At,  and  so  on.  The  phi's  and  theta's  are  constants,  is 
the  value  at  time  t  of  a  sequence  of  independent  random  impulses, 

9 

normally  distributed  with  variance  (white  noise),  a 
impulse  occurring  at  t-  At,  and  so  on. 


is  the 
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Comparing  Equations  4-1  and  4-4,  we  sec1  that  if  we  use  discrete 
Fourier  coefficients  to  estimate  a  weather  variable,  the  estimate  is 
deterministic;  it  is  a  function  only  of  time.  (The-  right-hand  side 
of  Equation  4-1  contains  only  constants  and  the  variable  i,  where  t  = 
i*t.)  If  we  use  an  ARMA  model  to  estimate  the  weather  variable,  the 
estimate  is  stochastic;  it  is  a  function  only  of  past  values  of  the 
variable  and  the  current  and  past  random  impulses  (the  estimate  is 
i ndependent  of  time) . 

A  third  analysis  technique  used  was  a  combined  deterministic 
plus  stochastic  model.  For  this  approach,  sinusoids  were  used  to 
account  for  the  deterministic  periodic  component  in  the  data,  and 
autoregressive  moving  average  terms  were  used  to  account  for  the  sto¬ 
chastic  component. 

We  will  show  that  the  results  of  the  application  of  these  three 
techniques,  taken  together,  allows  the  deterministic  periodic 
behavior  of  the  weather  data  to  be  identified. 

We  will  now  develop  a  few  of  the  ARMA  modeling  properties  that 
are  particularly  useful  for  weather  data  time  series  applications. 
More  complete  discussions  can  be  found  in  the  works  of  Kapoor  (ref 
28)  and  Box  and  Jenkins  (ref  26).  Pandit  (ref  29)  presents  a 
detailed  mathematical  and  historical  treatment  of  ARMA  modeling 
theory,  and  Kuo  (ref  30),  Franklin  and  Powell  (ref  31), 


and  Kline  (ref  32)  discuss  ARMA  modeling  in  the  context  of  digital 
control  theory. 

We  can  view  an  ARMA  model  as  a  discrete  linear  filter  (see  Fig¬ 
ure  4-1),  which  produces  a  time  series,  ,  as  output  when  subjected 
to  an  input  of  random  impulses  ^white  noise).  This  filter  has  vari¬ 
ous  forward  paths  and  -feedback  loops,  depending  on  the  order  of  t he 
ARMA  model.  This  visualization  helps  define  the  procedures  used  to 
find  an  adequate  ARMA  model  for  a  given  time  series.  We  invert  the 
ARMA  filter  as  part  of  the  model  identification  process.  By  viewing 
as  input,  we  can  construct  a  time  series,  a^ ,  for  any  candidate 
ARMA  model  (see  Figure  4-2).  This  procedure  is  made  obvious  by  rear¬ 
ranging  Equation  4-4: 


a.  =  Q,a.  ,  +  09at  „+'*'+  Oa.  „  +  X.  -  4* ,X+  , 
t  1  t-1  2  t-2  n  t-n  t  1  t-1 


-  <t>  x  -''‘-OX 

a  t-2  Vt-m 


[Eq  4-b] 


Given  that  we  can  calculate  a^,  there  are  two  steps  to  the  model 
identification  process: 


1.  For  a  model  of  a  given  order,  e.g.,  ARMA  (4,3),  we  must  find 
the  best  values  for  the  model  parameters  (the  phi's  and  theta's). 

This  is  accomplished  by  repeatedly  calculating  the  complete  time 
series  while  systematically  adjusting  model  parameters  until  the  sum 


of  squares  of  a^  is  minimized. 
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2.  Given  that  we  can  find  the  best  parameters  for  a  given  ARMA 
model,  we  further  require  that  the  a  time  series  produced  by  the 
"inverse"  ARMA  filter  be  random,  normally  distributed  impulses  (white 
noise).  Thus,  we  must  repeat  step  1  above  for  candidate  ARMA  models 
of  increasing  order  until  this  criterion  is  met. 

The  computer  software  system  entitled,  "University  of  Illinois 
Dynamic  Data  System  (UIDDS),"  developed  by  W.  A.  Kline  and  R.  L. 

DeVor  (ret  113),  provided  the  tool  needed  to  identify  adequate  ARMA 
models  for  the  weather  data  time  seriis  under  study.  The  first  step 
in  the  model  identification  process  is  largely  automatic  in  UIDDS 
through  the  use  of  optimisation  algorithms  which  usually  yield  the 
best  values  for  the  parameters  of  a  given  ARMA  model.  Part  2  of  the 
model  identification  process  is  not  automatic;  however,  selection  of 
r.n  adequate  moot1  from  ( <«nd  i  date  ARMA  models  is  made  possible  by  the 
wide  range  of  information  that  UIDDS  provides  for  each  candidate  ARMA 
model  . 


I  hr  [roceriun  tor  selecting  candidate  ARMA  models  recommended  by 
Piindil  (rtf  /'*)  w.  s  used  'o  analyze  the  weather  data  time  series  with 
UIDDS.  Inst,  parameters  for  an  ARMA  (P.l)  model  were  determined. 
Next.,  model  order  was  increased  according  to  ARMA  (iln.Kn-l),  n  =  2, 

?,  4...  until  an  adequate  model  was  obtained. 
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Model  adequacy  was  assessed  using  a  variety  of  statistics 
designed  to  establish  (1)  the  degree  to  which  the  a  lime  scries 
approximates  white  noise,  and  (2)  whether  increasing  the  model  order 
yields  statistically  "better"  models.  The  autocc rrel ati on  function 
and  the  Q-statistic  were  used  as  a  measure  of  the  randomness  of  the 
at  series.  The  autocorrel ation  function  is  a  series  which  measures 
the  correlation  of  a^  on  itself.  For  example,  the  value  of  the  auto¬ 
correlation  function  at.  lag  one,  p^,  is  simply  the  correlation 
between  values  of  the  residual  series  at  time  t  and  values  of  the 
series  at  t-1;  it  is  the  correlation  between  at  and  at  j.  Similarly, 
p^,  the  autocorrelation  at  lag  two,  is  the  correlation  between  a^  and 
at  2 *  and  so  on,  for  pk,  k  =  1,2, 3, 4, 5...  If  the  a  series  produced 
by  a  given  "inverted"  ARMA  model  approaches  white  noise,  then  the 
series  will  be  uncorrel ated ;  that  is,  values  of  the  autocorrelation 
function  at  all  lags  will  be  nearly  zero.  UIDDS  produces  plots  of 
the  autocorrel ation  function  which  were  inspected  for  each  candidate 
ARMA  model  to  determine  whether  the  values  of  the  function  were  all 
nearly  zero.  The  Q-statistic  is  a  cumulative  measure  of  autocorrela¬ 
tion;  it  is  simply  the  sum  of  squares  of  the  autocorrelation  function 
times  N,  the  number  of  points  in  the  time  series.  The  Q-statistic 
has  been  shown  to  vary  like  the  ciii-squared  statistic  and  can  be  com¬ 
pared  to  chi-squared  for  a  given  confidence  limit  and  degrees  of 
freedom  of  the  model . 
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The  F-test  and  Akaikt  Information  Criteria  (A1C)  were  used  as  a 
measure  of  model  "goodness."  (They  were  used  to  distinguish  between 
models  whose  autocorrelation  functions  arid  Q-stati sties  were  roughly 
the  same.)  The  F-test  is  a  comparative  test  between  two  models,  one 
of  which  contains  more  terms  than  the  other.  It  is  a  test  of  the 
hypothesis  that  the  extra  parameters  in  the  higher-order  model  are 
zero.  If  the  F-test  is  large,  the  parameters  are  probably  signifi¬ 
cant.  The  F-test  is  not  a  test  of  model  adequacy,  but  rather  an 
indicator  of  the  benefit  of  increasing  the  order  of  ARMA  models. 

AIC,  unlike  the  F-test,  is  a  measure  of  a  model's  absolute  merit.  It 
is  a  measure  of  the  model's  information  content  and  has  its  origins 
in  information  theory.  The  lower  the  AIC,  the  better  the  model;  how¬ 
ever,  a  model  may  bf  judged  to  be  adequate  even  if  other  models  have 
si ightl y  1 ower  AICs . 

Resul ts 

The  first  sets  of  weather  data  analyzed  were  hourly  measured 
dry-bulb  temperature  and  measured  solar  beam  radiation  for  Charles¬ 
ton,  SC,  for  the  year  1953.  The  first  40  days  of  the  data  were  used 
for  this  study;  for  most  of  this  first  modeling  effort,  every  other 
hour  of  the  data  set  was  removed  because  of  the  difficulty  in  apply¬ 
ing  ARMA  modeling  techniques  to  data  having  an  important  period  every 
24  sampling  intervals.  This  period  of  24  sampling  interval s--the 


diurnal  cycle  for  hourly  weather  data--was  too  long  to  be  success¬ 
fully  modeled  by  ARMA  methods. 

The  data  are  obviously  periodic  (If'  samples  per  24-hour  day). 
However,  because  solar  radiation  is  zero  at  night  (never  negative), 
beam  radiation  is  square-wave-like  when  compared  to  the  more 
sinusoidal  temperature  data.  One  way  of  approximating  a  square  wave' 
with  a  sinusoid  is  to  use  several  harmonics  of  the  fundamental  fre¬ 
quency.  Therefore,  we  anticipate  that  several  harmonics  of  the  diur¬ 
nal  frequency  may  be  needed  to  model  beam  radiation. 

We  begin  with  the  discrete  Fourier  spectrum  of  the  temperature 
and  the  beam  radiation  time  series  shown  in  Figures  4-3  and  4-4, 
respectively.  The  very  clear  peaks  at  40  cycles  per  sample,  or  one 
cycle  per  day,  and  at  the  first  few  harmonics  of  the  daily  cycle  are 
evident  for  both  temperature  and  radiation.  As  expected,  the  first 
harmonic  is  more  pronounced  for  the  beam  radiation  than  for  tempera¬ 
ture,  due  to  the  square  wave  nature  of  the  radiation  time  series. 
Higher  harmonics  are  present  in  both  time  series  but  are  compara¬ 
tively  insignificant. 

We  now  turn  to  results  of  the  ARMA  modeling  techniques  a ppl i ca¬ 
tion .  These  results  are  summarized  in  Tables  4-1  and  4-2,  which  show 
the  residual  sum  of  squares,  AIC,  the  Q  test,  the  various  F  tests, 
and  the  coefficient  of  determination,  R2,  for  various  ARMA  models  for 


both  time  series.  An  ARMA  (10,9)  model  proved  to  be  the  best  model 


Tabic  4-1 


Statistics 

for  Dry- 

-Bulb 

Temperature 

ARMA  Models 

F  Test  to 

RSS 

AIC 

Q/Q(95) 

Model  -  # 

R2 

1 

ARMA(4 ,3] 

1623 

60C 

202/38 

.89 

2 

ARMA(6 ,5) 

1507 

573 

170/38 

9-1 

.90 

2 

ARMA(8,7) 

1238 

487 

63/38 

25-2 

.92 

4 

ARMA( 10 ,9) 

1149 

459 

34/38 

9-3 

.92 

5 

ARMA( 12 , 11 ) 

1142 

465 

34/38 

.6-4 

.92 

6 

ARMA( 10,8) 

1179 

469 

45/38 

3-5 

.92 

7 

ARMA(8,6) 

1195 

468 

38/38 

-16-3 

.92 

Table  4-2 

Statistics  for  Solar  Beam  Radiation  ARMA  Models 


F  Test  to 


RSS 

AIC 

Q/Q( 95 ) 

Model  -  # 

R2 

1 

AP,MA(4 ,3) 

12.8E6 

4910 

135/38 

.59 

2 

ARMA(6,5) 

11.286 

4850 

53/38 

16-1 

.65 

? 

ARMA(8,7) 

10.3E6 

4819 

25/38 

9.8-2 

.67 

4 

ARMA{ 10,9) 

9.58E6 

4  793 

22/36 

8.5-3 

.70 

5 

ARKA( 12,11) 

9.72E6 

4808 

33/38 

-1.6-4 

.69 

6 

ARMA( 1C  .8) 

9.77E6 

4801 

23/38 

.69 

7 

ARMA (8, 6) 

10.1E6 

4806 

20/36 

.68 

3 

ARMA ( 6 ,4 ) 

10.4E6 

4816 

19/25 

.67 

for  both  temperature  and  solar  beam  radiation.  However,  in  both 
cases,  the  (Hi, 9)  model  had  a  small  0^,  and  the  (1(J,8)  model  was 
nearly  as  good,  particularly  for  solar  beam  radiation. 


Recall  that  we  have  visualized  an  ARMA  model  as  a  discrete 
filter  which  transforms  white  noise  into  the  requisite  time  series. 

By  applying  the  Z  transiorm  to  an  ARMA  (m,n)  model  of  Equation  4-4,  a 
discrete  transfer  function,  G(Z),  can  be  written  relating  X ( Z)  to 
a(Z),  where  X(Z)  is  the  Z  transform  of  X^ ,  and  a(Z)  is  the  Z 
transform  of  a^: 


C(Z) 


(1  -  0  7- 1  -  o  r2  -  e.J-3  -  •••  -  0  2"n ) 

_ _ 1 _ 2  o  _ n  ' 

(1  -  4'1Z'1  -  4>2Z-2  -  $  *  -  *3Z"ir') 


[Eq  4-6] 


Like  the  transfer  function  for  a  digital  filter  or  digital  con¬ 
troller,  this  transfer  function  has  real  end/or  complex  poles  and 
zeros.  Also,  since  at  is  assumed  to  be  white  noise,  the  theoretical 
spectrum  for  the  model  can  be  calculated  in  the  same  way  that  fre¬ 
quency  response  is  calculated  for  a  digital  filter  or  control ler--by 
finding  the  magnitude  of  G(Z)  for  various  frequencies  <o,  with  Z  set 
equal  to  e3"'  (.j  =  f-  1) .  Viewed  in  this  way,  an  autoregressive  moving 
average  model  can  be  analyzed  using  pole  placement  and  frequency 
response  techniques  common  to  digital  control  theory. 


For  example,  in  the  case  of  both  time  series  analyzed,  pairs  of 


complex  conjugate  poles  and  zeros  began  to  appear  near  the  unit 


circle  in  the  complex  Z-plain  tor  ARMA  models  of  order  as  low  as 
(4,3),  as  shown  in  the  sample  pole-zero  plots  of  Figure  4-b.  These 
pole  zero  pairs  were  of  nearly  equal  value.  Near,  but  not  exact, 
complex  conjugate  pole  zero  cancellation  near  the  unit  circle  is 
characteristic  of  periodic  time  series.  Fiowever,  for  both  series, 
two  poles  remained  well  within  the  unit  circle  as  the  order  of  the 
models  was  increased.  One  zero  tended  to  appear  near  the  origin  in 
the  higher-order  ARMA  (2n,  2n-l)  models.  The  conclusion  which  can  be 
drawn  from  this  behavior  is  that  the  time  series  are  characteri zed  by 
periodic  behavior  (indicated  by  the  pole  zero  pairs  near  the  unit 
circle)  and  by  stochastic  AR(2)  behavior  (indicated  by  the  remaining 
poles  which  are  non-zero  but  well  inside  the  unit  circle). 

We  can  also  compare  the  spectrum  of  the  autoregressive  moving 
average  model  for  temperature  and  solar  radiation  to  the  Fourier 
spectrum  computed  previously.  Figures  4-6  and  4-7  show  the  spectrum 
of  the  ARMA  (10,9)  models  for  temperature  and  radiation,  respec¬ 
tively.  We  immediately  see  one  difficulty  which  arises  when  strongly 
periodic  functions  are  modeled  using  stochastic  ARMA  methods.  The 
spectral  peaks  for  the  ARMA  transfer  functions  are  at  the  same  fre¬ 
quencies  as  those  of  the  discrete  Fourier  transform  (see  Figures  4-3 
and  4-4).  However,  the  ARMA  spectrum  is  a  continuous  spectrum,  and 
we  must  integrate  the  square  of  the  spectrum  shown  in  Figures  4-6  or 
4-7  from  to-m/N  to  ui+ti/N  in  order  to  compare  the  "power"  at  frequency 
ix i  with  the  "power"  represented  by  the  square  of  the  discrete  line 
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Figure  4-7.  Model  spectrum--beam  solar  radiation.  Charleston,  SC 


spectrum  at  m.  While  we  can  carry  out  this  integration,  it  is  tedi¬ 
ous.  There  are  also  other  mathematical  difficulties.  If,  as  physi¬ 
cal  reasoning  and  the  discrete  Fourier  spectrum  suggest,  the  weather 
data  contain  deterministic  sinusoidal  components,  then  the  continuous 
spectrum  of  a  climate  model  containing  these  sinusoids  would  be 
discontinuous  at  the  frequencies  corresponding  to  the  principal 
sinusoids.  In  order  for  the  ARMA  spectrum  to  have  such  discontinui¬ 
ties,  pairs  of  conjugate  poles  would  have  to  appear  on  the  unit  cir¬ 
cle,  not  just  near  it.  Such  a  model  would  not  be  stable  and  would  be 
of  little  utility.  If  we  are  careful  to  choose  our  sampling  interval 
so  that  the  periods  of  the  sinusoids  in  the  data  are  exact  integer 
multiples  of  the  s&mpling  interval,  then  the  discrete  Fourier 
transform  provides  the  best  estimate  of  the  amplitude  of  the 
sinusoids.  The  sharp  spikes  at  exactly  the  diurnal  frequency  and  i 
first  harmonic  in  the  ARMA  spectrum  confirm  that  periodici ty  exists  at 
exactly  these  frequencies  which  are  integer  multiples  of  the  sampling 
interval . 

We  must  be  satisfied  that  ARMA  modeling  alone  confirms  the  need 
for  combined  deterministic  plus  stochastic  models.* 

We  deal  briefly  now  with  the  reason  for  removing  every  other 
hourly  sample  from  this  original  weather  data  set.  If  we  consider  a 
simple  sinusoid  with  24  samples  per  period,  we  see  that  much  of  the 
★ 


We  use  the  phrase  "deterministic  plus  stochastic  model"  to  identify 
a  single  model  with  deterministic  and  stochastic  components. 
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sine  wave  is  like  a  straight  line  from  the  microscopic  view  of  one 
sample  looking  at  its  nearest  four  or  five  neighbors  (from  the  per¬ 
spective  of  an  ARMA  model).  Thus,  we  can  estimate  the  value  of  the 
time  series  at  time  t  with  reasonable  accuracy  by  looking  only  at  a 
few  past  values  of  the  time  series  and  one  or  two  values  of  the  past 
a^  time  series.  For  such  a  finely  sampled  time  series,  the  residual 
sum  of  squares  can  be  reduced  substantial ly  with  a  model  that  does 
not  portray  the  periodic  behavior  of  the  data  at  all. 

The  results  of  the  spectral  analysis  and  ARMA  modeling  both  sug¬ 
gest  the  need  for  a  deterministic  (i.e.,  time-dependent  sinusoids) 
plus  stochastic  (i.e.,  low-order  ARMA)  model.  The  U1DDS  software 
system  also  allows  the  parameters  of  combined  deterministic  plus  sto¬ 
chastic  models  to  be  estimated  using  nonlinear  regression. 

To  begin  this  modeling  step,  the  form  of  a  combined  model  and 
the  initial  estimates  of  the  model  parameters  must  be  established  by 
first  finding  an  adequate  deterministic  only  model  and  then  finding 
an  adequate  ARMA  model  for  the  residuals  left  after  the  deterministic 
behavior  is  removed.  Thus,  for  both  the  temperature  and  beam  radia¬ 
tion  time  series,  sinusoidal  deterministic  models  were  tested  by 
first  including  just  the  diurnal  frequency,  then  the  first  harmonic, 
then  another  harmonic,  and  so  on,  until  there  was  no  substantial 
improvement  in  the  residual  sum  of  squares  between  successive  models. 
Tables  4-3  and  4-4  show  the  results  of  these  tests.  In  both  cases. 


Table  4-3 


Deterministic  Model  Statistics,  Dry-Bulb  Temperature,  Charleston,  SC 

Standard 


Modt  1 

Periods 

Amp  I  itudes 

KSS 

R2 

Deviation 

i 

13 

4.U/ 

11.4313 

.Zb 

4.88 

> 

13 

4.0? 

1 1 .00t 3 

.29 

4.79 

b 

1.33 

j 

13 

4.0/ 

10.9813 

.29 

4.78 

0  1.33 

3  .3  a 

Standard  deviation  ot  tin-  data  -  S.66 


Table  4-4 

Deterministic  Model  Statistics,  Solar  Beam  Radiation,  Charleston,  SC 


Standard 


Mudt*l 

Periods 

Amp! itudes 

RSS 

R2 

Deviation 

i 

12 

219.0 

20.1 01  (. 

.  30 

204 

2 

12 

219.0 

Hi. SHI  0 

.48 

IMS 

6 

120.0 

3 

12 

219.0 

16.4710 

.48 

18b 

6 

120.6 

2.4 

20.8 

Standard  deviation  ot  the  data  *  2 S7 


the  fundamental  diurnal  frequency  and  its  first  harmonic  were  needed, 
but  use  of  additional  terms  did  not  significantly  improve  the  deter¬ 
ministic  models.  Notice  that  the  addition  of  the  first  harmonic  pro¬ 
duced  a  much  larger  change  in  the  residual  sum  of  squares  for  the 
solar  beam  radiation  series  than  for  the  temperature  series;  this 
confirms  our  earlier  expectation  that  the  first  harmonic  would  play  a 
more  important  part  in  the  characterization  of  the  square-wave-like 
solar  data. 

Next,  adequate  ARMA  models  were  developed  for  the  residuals  from 
two-term  deterministic  models.  For  each  of  the  two  series,  an  AR(2) 
model  was  found  to  be  adequate. 

Finally,  the  model  parameters  were  estimated  for  combined 
sinusoidal/AR(2)  models  for  both  series.  Tables  4-6  and  4-6  show 
the  final  results.  Note  that  even  though  two  sinusoidal  terms  were 
found  to  be  significant  in  the  deterministic  phase  of  the  model 
development,  in  the  combined  phase,  three  sinusoids  plus  AR ( 2 )  models 
proved  best  for  both  weather  data  sequences,  although  the  third 
sinusoid  has  a  very  small  amplitude.  The  residual  sum  of  squares  and 
the  A I C  for  these  best  combined  models  is  lower  than  their  counter¬ 
parts  for  the  best  ARMA  models. 

We  now  compare  the  results  of  the  three  analysis  techniques 


applied  to  the  two  Charleston  weather  data  time  series. 
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Table  4-5 

Deterministic  Plus  Stochastic  Models  for  Every  Other  Hour  Dry-Bulb 
Temperature,  Charleston,  SC,  1953 


Model 

Peri od 

Ampl itude 

Residual 
Sum  of 
Squares 

A1C 

R2 

Standard 

Deviation 

Standard 
Deviation 
of  the  Data 

1 

12 

4.07 

1229 

4S5 

.92 

1 .60 

5.66 

6 

1.33 

PHI (1 )  =  1.183 

PHI (2)  =  - .269 

2 

12 

4.06 

1128 

414 

.93 

1  .S3 

5.66 

6 

1.33 

3 

.33 

PH 1 ( 1 )  =  1.230 
PHI (2)  =  -.303 


Table  4-6 

Deterministic  Plus  Stochastic  Models  for  Every  Other  Hour  Solar 

Beam  Radiation 


Residual  Standard 

Sura  of  Standard  Deviation 


Model 

Pt  r  i  od 

Ampl i tude 

Squares 

A1C 

R2 

Deviation 

of  the  Data 

1 

12 

6 

PHI(l)  =  .746 

PH  1 (2)  -  -.1/8 

220.1 

120.4 

9.612  x  106 

4758 

.70 

141 

257 

2 

12 

220.4 

9.261  x  106 

4740 

.71 

139 

257 

6  120.2 

2.4  20. H 

I’M  1(1)  -  .// 4 
PHI  (2)  -  - .2(H) 


Both  the  discrete  Fourier  transform  spectrum  and  the  ARMA  model¬ 
ing  procedure  allow  identification  of  deterministic  periodic  behavior 


lib 


in  the  weather  data  time  series  studied.  However,  even  t hough  ade¬ 
quate  ARK A  models  can  be  found  for  the  strongly  periodic  time  series, 
the  ARMA  model  spectra  are  not  adequate  indicators  of  spectral  peak 
heights  in  the  data;  consequently,  they  do  not  provide  good  estimates 
of  the  amplitudes  of  the  principal  sinusoids  in  the  data.  Notice, 
however,  that  the  amplitudes  of  the  sinusoidal  terms  of  the  combined 
deterministic  plus  stochastic  models  are  very  nearly  equal  to  the 
height  of  the  spectral  peaks  of  the  discrete  Fourier  transform  (see 
Figures  4-1  and  4-2).  For  the  time  series  studied,  and  for  other 
deterministic  periodic  time  series,  the  discrete  Fourier  transform 
can  provide  good  amplitude  estimates. 

Generally,  the  combined  model  is  more  satisfactory  than  ARMA- 
only  models  for  characterizing  weather  data  time  series,  since  the 
deterministic  periodic  behavior  of  the  data  is  separated  from  the 
stochastic  autoregressive  behavior. 

The  development  of  pure  ARMA  models  is  unnecessary  for  these 
weather  data  time  series,  since  their  discrete  Fourier  spectral  peaks 
are  extremely  pronounced.  A  combined  model  can  be  developed  directly 
from  amplitude  and  period  information  obtained  from  the  Fourier  spec¬ 
trum  and  from  ARMA  modeling  of  the  residuals.  We  can  carry  out  this 
procedure  on  a  sequence  of  hourly  data  rather  than  on  data  from  every 
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('Mur  hour-,  since  t.hr  problems  associated  with  ARMA-only  modeling  ol 
linrly  sampled  dot  a  sequence  described  earlier  will  no  longer  be 

r : ■  1 ev, n  t  . 


Wc  must  ensure  that  we  correctly  interpret  the  discrete  Fourier 
spectrum  and  the  spectra  of  the  ARMA  models  at  points  away  from  the 
spikes  in  the  spectrum.  The  small  "hill"  beginning  at  zero  frequency 
and  tailing  off  rapidly  in  the  temperature  data  spectrum  is  charac¬ 
teristic  of  an  AR(2)  model;  it  does  not  represent  the  presence  of 
low-frequency  components  in  the  weather  data  series.  From  the  fre¬ 
quency  perspective,  only  the  diurnal  frequency  and  its  harmonics 
exist,  in  the  weather  data  time  series.  These  are  the  only  frequen¬ 
cies  which  are  physically  justifiable.  (Longer  periods  may  exist  in 
longer  data  records.) 

To  verify  that  our  conclusion  from  the  Charleston  data  analysis 
can  be  generalized,  Fourier  spectra  were  plotted  for  a  number  of  sam¬ 
ples  from  different  sites  during  different  times  of  the  year.  These 
spectra,  shown  in  Figures  4-8  through  4-13,  are  for  records  of  hourly 
samples  (/?(!  hours  per  sample)  rather  than  samples  taken  every  other 
hour.  In  addition,  deterministic  plus  stochastic  models  were 
developed,  using  the  diurnal  and  the  first  harmonic  of  the  diurnal  as 
the  frequency  components  of  the  deterministic  model.  In  each  case, 
the  autoregressive  component  was  AR(2).  Notice  that  the  AR(2) 
behavior  persists,  even  though  we  have  changed  the  sampling  interval. 


Figure  4-10.  Spectrum  for  temperature,  Madison,  WI 
(Zb  September  through  25  October  1955). 
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Figure  4-11.  Spectrum  for  beam  solar  radiation,  Madison,  WI 
(25  September  through  25  October  1960). 


table  4-/  shews  L hr  coefficients  of  deterministic  plus  stochastic 
models  for  Madison,  WI,  Fort  Worth,  TX,  and  Santa  Maria,  CA. 

having  characterized  the  short-term  behavior  of  weather  data,  we 
must  determine  tin  extent  of  the  annual  weather  cycle  (which  we  know 
i x i s" s  on  physical  grounds)  and  determine  if  there  is  any  other 
periodic  behavior  in  the  data  between  the  diurnal  and  the  annual 
cycles.  We  can  also  look  for  stochastic  behavior'  in  longer  aata 
records.  To  examine  long-term  weather  data  trends,  two  time  series 
were  constructed  for  both  temperature  and  beam  radiation;  one  con¬ 
sisted  of  the  daily  averages  of  the  hourly  data  for  a  Ill-year  period 
(3f>50  points),  the  other  consisted  of  b-day  averages  of  the  hourly 
dal  a  lor  a  li'-year  period  (  730  points).  Data  from  Charleston,  SC, 
Madison,  WI,  and  Santa  Maria,  CA,  were  used  for  this  part  of  the 
s' udy. 

As  before,  we  first  examine  the  Fourier  spectra  of  each  time 
series.  The  time  scries  of  daily  averages  is  used  to  ensure  that  any 
cyclical  behavior  with  a  period  between  2  and  10  days  will  show  up 
clearly.  Any  reguK»  behavior  with  periods  longer  than  10  days  will 
b<  revealed  by  the  spectra  of  the  time  series  of  b-day  averages; 
these  lime  series  arc  also  short  enough  to  be  analyzed  using  the 
deterministic  plus  stochastic  modeling  techniques. 


Table  4-7 

Deterministic  Plus  Stochastic  Models  for  Various  Weather 

Data  Time  Series 


Weather 

Stochast  ic 

Residual  Sum 

Standard 

Standard 

Local  ion 

Data 

Period 

Arnpl  Itude 

Parameters 

of  Squares 

R2 

Dev  la? Ion 

Devlat (on  of  Dai  a 

fort  Worth,  IX 

temperature 

24 

4.04 

PH  1(1) 

-  1.363 

743 

.98 

1.01 

8.41 

1  feb  -  2  Mar  55 

12 

1.05 

PHI (2) 

-  -.3/7 

Solar  Be dm 

2 A 

296 

PHI (1) 

*  .92/ 

'  1 .0/0  *  10? 

.86 

122 

324  1  1 

Radiat ion 

1? 

114 

PHI (2) 

•  -.ORB 

| 

Madison,  Wl 

Temperature 

24 

4.49 

PHI(n 

*  1.178 

1004 

.97 

1.18 

6.62  !’ 

25  S*-|i  2b  (Jet  6U 

12 

1.49 

PHI (? 1 

-  -.206 

j 

Solar  Beam 

24 

21 9  * 

PHI { 1 ) 

-  .863 

.984  x  10> 

.84 

117 

288  i  | 

Radiation 

12 

124 

PHI (2) 

-  -.096 

Santa  Mar  Id.  (.A 

Temperature 

24 

4.04 

PHI  U  1 

-  1.0/5 

385 

.95 

.73 

3.23 

29  Jun  28  Jut  hi 

1, 

MB 

PHI  (2 ) 

-  -.283 

Solar  Beam 

24 

458 

PHI (1 1 

-  .94  5 

.725  x  10/ 

.93 

100 

373 

Radial  ion 

1? 

144 

PH  1  (?) 

-  -.291 

1 
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Figure  4-14  shows  the  spectrum  for  the  daily  average  temperature 
time  series  for  Madison,  WI.  There  is  quite  clearly  only  one  spec- 
'ral  peek,  which  occurs  at  1C  cycles  per  sample  or  one  cycle  per 
year,  everywhere  else,  the  spectrum  resembles  a  "white  noise"  spec¬ 
trum  (at  least  by  comparison  to  the  annual  cycle  peak).  We  cannot 
yet  conclude  that  all  that  will  remain  after  the  annual  cycle  is 
removed  will  be  random  noise,  but  Figure  4-14  clearly  allows  us  to 
eliminate  short-term  (2-  to  10-day)  periodicity  from  further  con¬ 
sideration.  This  was  confirmed  by  the  spectre  for  the  other  five 
time  series  of  daily  averages  (not  shown)  which  had  exactly  the  same 
spectral  claracter. 

We  now  examine  the  spectra  for  the  time  series  of  6-day  averages 
(Figures  4-16  through  4-20).  Again,  the  annual  cycle  (1C  cycles  per 
sample)  is  the  only  obvious  characteristic  of  these  data.  The  first 
harmonic  of  the  annual  cycle  appears  only  for  beam  radiation  for 
Charleston,  SC;  in  this  case,  both  the  annual  cycle  and  its  first 
harmonic  are  relatively  small  when  compared  to  the  mean.  On  the 
basis  of  these  spectra,  we  conclude  that  a  combined  deterministic 
plus  stochastic  model  for  the  5-day  average  time  series  should 
include  only  the  annual  sinusoid  as  its  deterministic  component. 

To  determine  whether  stochastic  behavior  (other  than  random 
variance)  existed  in  the  weather  data  time  series,  ARMA  modeling 
techniques  were  applied  to  the  residuals  from  a  deterministic  model 
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Figure  4-17.  Spectrum  for  5-day  average  temperature,  Charleston, 
SC  (1953-1962). 


CYCLES  PER  SAMPLE 


Figure  4-18.  Spectrum  for  5-day  average  beam  solar  radiation 
Charleston,  SC  (1953-1962). 


for  each  time  series  (the  deterministic  model  included  only  the 
annual  sinusoid).  An  AR(1)  model  was  found  to  be  adequate  in  each 
case.  Finally,  combined  deterministic  plus  stochastic  models  were 
found  for  each  time  series. 

Table  4-8  shows  the  results  for  the  deterministic  only  models 
and  for  the  deterministic  plus  stochastic  models.  We  see  that  for 
climates  like  Madison  and  Charleston,  which  have  pronounced  annual 
temperature  variations,  the  annual  sinusoid  accounts  for  a  large 
fraction  of  the  variance  in  the  temperature  data.  For  Santa  Maria, 
annual  temperature  swings  are  small,  and  the  annual  temperature  cycle 
is  less  important  than  daily  temperature  cycles  (see  Tables  4-6  and 

4- 7  for  compari son) . 

Comparing  the  annual  cycles  for  solar  beam  radiation  with  the 
daily  cycles  of  Tables  4-6  and  4-7,  we  see,  as  expected,  that  annual 
variations  are  much  less  pronounced  than  daily  variations.  This  is 
particularly  true  for  Charleston  with  its  southern  latitude. 

We  can  assess  the  importance  of  the  AR  components  in  the  models 
shown  in  Table  4-8  by  comparing  the  standard  deviation  of  the  deter¬ 
ministic  only  models  with  its  counterpart  for  the  deterministic  plus 
stochastic  models.  At  each  site,  we  see  that  the  AR(1)  component 
produces  only  a  nominal  reduction  in  the  standard  deviation  for  the 

5- day  average  temperature  data  (less  than  .2°C)  and  almost  no  reduc¬ 
tion  for  the  5-day  average  solar  data.  The  very  small  improvement  in 
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Table  4-8 

Deterministic  and  Deterministic  Plus  Stochastic  Models  for 
5-Day  Averaged  Weather  Data  Time  Series 


f  requenry 


l  (k  at  ion 

Oat  a  {( ye  Its/ Year) 

Amp I i t  udr 

PHI { 1  ) 

Squar es 

lit  v  i at  l 
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3.04 
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bj .  UU 
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Temperature  1 
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24bb 

.bb 

l.U 
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.(.4 

1.64 

Solar  Beam  1 

Radial  ion 
(Del ermini st i(  ) 

66.1K. 
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Solar  Re, mi  1 
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Plus  Stochastic) 

6b.  9b 

.2/2 

2.21b  x  1111’ 

,4  a 

bs.nP 

the  models  due  to  the  AR{1)  component  does  not  justify  the  additional 
model  complexity  that  including  this  term  would  cause. 


Having  analyzed  both  the  short-term  (hourly)  and  long-term  (b- 
day  averages)  character  of  weather  data  time  series,  we  can  now  write 
the  general  form  for  a  deterministic  plus  stochastic  model  with  a  1- 
hour  time  step  for  both  temperature  and  beam  solar  radiation  (models 
for  both  temperature  and  radiation  have  the  same  form) : 

xt  =  ao  +  2(ax  cos(2iTt/8760)  +  bL  sin(2  t/8760) ) 

+  2(a365  cos(2ir365t/8760)  +  b36b  si n(21T36bt/8760) ) 

+  2(ay30  cos(2tt730/8760)  +  b?30  si  n ( 2 -rT730t /8760 ) ) 

+  (other  frequency  components;  see  discussion) 


+  *lXt-l  +  Vt-2  +  at 


[Eq  4-7] 


where:  Xt  is  the  value  of  the  weather  data  variable  at  time  t 

a  is  the  mean  of  the  weather  data  variable 
o 

a^  and  b^  are  the  annual  cycle  Fourier  coefficients 


a 36 5  and  b365 


are  the  daily  cycle  Fourier  coefficients 


a73Q  and  b73g  are  the  Fourier  coefficients  of  the  first 


harmonic  of  the  daily  cycle 


and  ^  are  constant  autoregressive  model  parameters 
i  and  2  are  the  values  of  the  weather  data  variable  at 
t-1  and  t-2,  respectively 


is  the  random  impulse  at  time  t. 
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Notice  that  we  have  used  the  Fourier  coefficients  as  the  determinis¬ 
tic  model  components  rather  than  coefficients  derived  from  nonlinear 
regressions.  This  is  because  the  amplitudes  of  the  sinusoids,  which 
result  from  the  nonlinear  regression  used  to  find  deterministic  plus 
stochastic  models,  are  almost  exactly  equal  to  the  corresponding 
Fourier  spectrum  peak  heights  for  all  the  cases  we  have  studied. 

Thus,  we  may  as  well  use  the  more  efficient  fast  Fourier  transform 
technique  to  find  the  deterministic  model  parameters. 

The  "other  frequency  components"  of  Equation  4-7  allow  for  the 
inclusion  of  other  harmonics  if  they  are  determined  to  be  significant 
for  a  particular  site  or  climate  variable  and  for  the  inclusion  of 
any  significant  side  bands.  Side  bands  were  shown  to  be  important  in 
the  work  described  in  ref  2  and  by  the  Fourier  transforms  of  1  year 
of  hourly  data  described  below  and  in  Chapter  5. 

Side  bands  are  spectral  peaks  which  occur  on  either  side  of  a 
fundamental  frequency  or  its  harmonics.  For  example,  spectral  peaks 
at  364  and  366  cycles  per  year  are  side  bands  of  the  diurnal  cycle 
(365  cycles  per  year  is  the  diurnal  cycle).  Side  bands  are  properly 
interpreted  as  indicating  the  presence  of  amplitude  modulation.  For 
example,  side  bands  at  364  and  366  indicate  that  the  amplitude  of  the 
diurnal  cycle  varies  sinusoidally  about  its  mean  value  (as  indicated 
by  the  height  of  the  spectral  peak  at  365  cycles  per  year).  The  fact 
that  the  side  bands  occur  at  365+1  indicates  that  the  amplitude  of 
the  diurnal  cycle  modulates  at  1  cycle  per  year. 

i 

_ i  i 
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Side  bands  around  the  diurnal  cycle  and  its  harmonics  are  most 
pronounced  for  solar  radiation  in  northern  climates  (e.g.,  Madison). 
This  simply  indicates  that  there  is  substantial  seasonal  variation  in 
the  amplitude  of  the  diurnal  variation  in  solar  radiation,  a  fact 
well  known  from  astronomy. 

Equation  4-7  represents  the  culmination  of  the  analysis  which 
has  accomplished  the  following:  separation  of  the  deterministic 
periodic  behavior  from  the  stochastic  behavior  in  the  weather  data 
time  series;  demonstration  that  there  are  only  a  few  statistically 
significant  frequencies  in  the  climate  data  (the  annual  cycle,  the 
diurnal  cycle,  and  the  first  harmonic  of  the  diurnal  are  the  most 
important;  other  harmonics  and  side  bands  may  be  important,  depending 
on  the  site  and  the  climate  variable);  and  discovery  that  the  weather 
data  time  series  has  a  significant  stochastic  component. 

As  a  final  step  in  the  weather  data  analysis,  we  need  to  deter¬ 
mine  sample  statistics  for  an  hourly  model  which  includes  the  mean, 
the  annual  cycle,  the  diurnal  cycle,  and  the  diurnal  cycle's  first 
harmonic,  but  does  not  include  the  stochastic  component.  To  achieve 
this,  fast  Fourier  transform  coefficients  were  determined  for  1  year 
of  hourly  data  for  1955  from  Fort  Worth,  Santa  Maria,  Charleston,  and 
Madison.  Hourly  temperature  and  solar  beam  radiation  were  modeled 


I 
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tor  rath  site  tor  ..  1-year  period,  using  only  the  mean,  the  annual 
cycle,  the  diurnal  cycle,  and  the  first  harmonic  ot  Un  diurnal 
cycle;  results  were  then  compared  to  the  actual  data  (see  Table  4-9). 
We  can  see  from  the  table  that  the  larger  deterministic  components  in 
the  weather  data  account  for  a  substantial  portion  of  the  variance  in 
the  data--more  of  the  variance  than  is  indicated  in  the  results  shown 
in  Table  4-3 . 

Corel  us  ion 

As  demonstrated  in  Chapter  3,  the  use  of  the  sinusoidal  com¬ 
ponents  of  the  weather  data  model  without  the  use  of  the  stochastic 
component  greatly  simplifies  the  calculation  of  loads  for  building 
energy  analysis.  The  analysis  presented  in  this  chapter  shows  that 
true  periodicity  in  the  data  exists  only  at  the  annual  cycle  and  at 
the  diurnal  cycle  plus  its  harmonics.  We  have  also  shown  that  the 
fast.  Fourier  transform  is  an  adequate  procedure  for  determining  the 
amplitude  of  the  periodic  components.  Using  the  fast  Fourier 
transform,  we  can  also  analyze  much  longer  periods  of  hourly  data 
than  we  have  shown  in  the  preceding  analysis  or  than  could  be 
analyzed  by  other  methods,  such  as  nonlinear  regression. 

These  weather  data  models  also  have  other  significant  applica¬ 
tions.  One  of  the  most  obvious  is  the  potential  use  of  deterministic 
plus  stochastic  models  in  lieu  of  weather  data  tapes  as  input  to 
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Table  4-9 

Deterministic  Sinusoidal  Models  for  Hourly  Weather 
Data  for  1  Year 


Loeat ion 

Weather 

Data 

R2 

Standard 

Dev iat ion 

Standard 
Deviation 
of  the  Data 

Fort  Worth,  IX 

Temperature 

.73 

5.61 

9.86 

Solar  Beam  Radiation 

.60 
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Charleston,  SC 

Temperature 

.72 

4.40 

8.32 

Solar  Beam  Radiation 

.60 

190 

269 

Madison,  W1 

Temperature 

.83 

5.36 

13.04 

Solar  Beam  Radiation 

.47 

203 

280 

Santa  Maria,  CA 

Temperature 

.59 

3.15 

5.94 

Sol ar  Beam  Radiation 

.66 

194 

335 

3 


'i 

i 
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tnergv  analysis  programs.  Equations  like  Equation  4-7  could  simply 
be  included  in  the  energy  analysis  programs  themselves,  and  the  model 
coefficients  would  be  part  of  the  input  to  the  program.  The  hourly 
value  for  a  could  be  produced  by  a  seeded  random-number  generator 
with  the  variance  of  a  supplied  as  input.  This  approach  could  com¬ 
pletely  eliminate  the  need  to  routinely  access  weather  data  tapes  or 
files  (simultaneously  eliminating  controversy  over  data  formats  and 
weather  data  file/energy  analysis  program  interfaces).  The  coeffi¬ 
cients  for  a  combined  deterministic  plus  stochastic  model  need  to  be 
generated  and  tabulated  only  once  for  any  given  site,  and  they  could 
be  based  on  very  long  periods  of  hourly  recorded  data.  The  genera¬ 
tion  of  the  coefficients  would  be  quite  straightforward.  The 
sinusoidal  terms  can  be  derived  using  fast  Fourier  transform  tech¬ 
niques  (available  on  most  computer  systems)  which  are  very  efficient 
and  can  accommodate  very  large  data  sets  (several  years  or  even  tens 
of  years  of  hourly  data).  The  stochastic  model  component  has  only  AR 
behavior,  which  means  that  the  model  parameter  can  be  determined 
using  1 i near  regression,  again  allowing  the  analysis  of  very  large 
data  sets. 

The  weather  data  models  which  would  result  from  the  above  data 
analysis  would  be  truly  "typical”  in  a  statistical  sense,  thereby 
perhaps  laying  to  rest  much  of  t.he  controversy  over  what  weather  data 
to  use  for  energy  analysis  purposes. 


Another  possible  application  of  the  deterministic  plus  stochas 


tic  climate  models  is  short-term  weather  forecasts  for  use  in 
microprocessor-based  building  energy  system  controls.  The  control  of 
chilled  water  storage,  for  example,  might  be  based  on  an  algorithm 
which  included  forecasts  of  outdoor  temperatures  and/or  solar  radia¬ 
tion  for  the  next  several  hours  or  days.  The  forecasts  would  be 
updated  with  measured  data  each  hour.  The  availability  of  a  good 
short-term  forecast  (perhaps  10  to  20  hours  ahead)  could  be  useful  in 
controlling  a  variety  of  KVAC  system  parameters,  such  as  morning 
startup  time,  chilled  and  hot  water  temperatures,  length  of  duty 
cycling  used  to  avoid  peak  loads,  and  staging  or  load  sharing  in 
large  boiler  or  chiller  plants.  The  use  of  climate  forecasts,  com¬ 
bined  with  the  calculated  or  measured  dynamics  of  the  building  and 
its  energy  systems  (room  and  system  "response  factors"),  would  allow 
the  control  system  to  anticipate  heating  and  cooling  requirements. 

The  potential  for  reducing  building  energy  consumption  by  applying 
control  systems  which  use  climate  forecasts  remains  to  be  explored. 
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5  LFFEC 13  OF  DL TERM  1 N I  ST  I C  CLIMATE  MODELS 
ON  ENERGY  CALCULATIONS 

Problem  Definition 

In  Chapter  3,  we  showed  that  calculating  heat  flow  through 
multi-layered  slabs  is  greatly  simplified  if  the  temperature  boundary 
condition  on  one  side  of  the  slab  is  assumed  to  be  sinusoidal  rather 
than  a  discrete  temperature  time  series.  In  particular,  calculating 
the  frequency-dependent  U-val ue  of  a  multi-layered  slab  is  much 
easier  than  calculating  conduction  transfer  functions.  We  expected 
to  be  able  to  use  frequency  response  methods  as  part  of  a  simplified 
p, ocedure  for  calculating  hourly  heating  and  cooling  loads  in  build¬ 
ings. 


We  then  sought  to  determine  how  many  terms  of  a  Fourier  series 
would  be  required  to  provide  an  "adequate"  approximation  of  an  actual 
weather  data  time  series.  The  results  of  the  application  of  time 
series  analysis  techniques  to  weather  data,  summarized  in  Chapter  4, 
led  us  to  conclude  that  truly  sinusoidal  behavior  is  exhibited  only 
at  very  few  frequencies  (one  cycle  per  year  and  one  cycle  per  day 
plus  harmonics  of  one  cycle  per  day).  The  rest  of  the  variation  in 
weather  data  is  stochastic  in  nature.  We  could  account  for  more 
variance  in  a  weather  data  time  series  by  using  greater  numbers  of 
Fourier  coefficients  to  approximate  the  series;  however,  using 
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coefficients  whose  frequencies  did  not  correspond  to  the  annual  cycle 
or  the  diurnal  cycle  and  its  harmonics  or  to  their  side  bands  would 
be  an  artificial  means  of  accounting  for  the  portion  of  climate  vari¬ 
ation  that  is  really  stochastic  rather  than  deterministic. 

In  this  chapter,  we  examine  results  of  a  study  to  determine  how 
much  error  will  be  introduced  in  calculating  expected  energy  use  in  a 
building  if  we  use  only  the  mean,  the  annual  cycle,  and  the  diurnal 
cycle  and  its  first  harmonic  to  characterize  weather  variation. 

The  procedure  used  to  carry  out  this  investigation  consisted  of 
the  following  steps: 

1.  For  each  of  four  representative  locations,  a  fast  Fourier 
transform  was  performed  on  1  year  of  hourly  weather  data  for  each 
climate  variable  affecting  building  energy  use  calculations. 


2.  The  mean  of  each  climate  variable  and  the  Fourier  coeffi¬ 
cients  corresponding  to  the  annual  cycle,  the  diurnal  cycle,  and  the 
first  harmonic  of  the  diurnal  cycle  were  used  to  synthesize  1  year  of 
hourly  weather  data. 

3.  The  energy  use  of  typical  buildings  served  by  several  dif¬ 
ferent  heating  and  cooling  systems  was  calculated  using  both  actual 
and  synthesized  hourly  weather  data.  The  Building  Loads  Analysis  and 
System  Thermodynamics  (BLAST)  program,  a  detailed  energy  analysis 
computer  program,  was  used  to  calculate  expected  energy  use. 


4.  Results  of  Step  3  were  compared  and  analyzed. 

Weather  Data 

The  same  weather  data  tapes  used  for  the  analysis  described  in 
Chapter  4  were  used  for  this  part  of  the  study.  Recall  that  the  four 
sites  selected  were  Charleston,  SC,  a  hot  humid  climate;  Madison,  WI, 
a  cold  northern  climate;  Fort  Worth,  TX,  a  hot  dry  climate;  and  Santa 
Maria,  CA,  a  marine  climate.  The  year  1955  was  selected  from  the 
available  data  records  as  being  reasonably  representative. 

The  BLAST  program  used  the  following  climate  variables  to  calcu¬ 
late  expected  building  energy: 

1.  Outdoor  dry-bulb  temperature 

?.  Outdoor  wet-bulb  temperature 

2.  Barometric  pressure 

4.  Humidity  ratio  (derived  from  dry-bulb,  wet-bulb,  and 
barometric  pressure  data) 

5.  Wind  speed 

6.  Beam  radiation 


7.  Diffuse  radiation 
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8.  Ground-reflected  radiation  (calculated  from  beam  and  diffuse 
radiation,  ground  reflectivity,  and  a  snow  cover  indicator). 

A  fast  Fourier  transform  was  performed  on  1  year  of  hourly  data 
for  each  climate  variable.  Figure  5-1  shows  sample  results  of  the 
Fourier  transform  for  beam  radiation  at  Santa  Monica,  CA.  Cnee  the 
fast  Fourier  transform  was  completed,  coefficients  representing  the 
mean,  the  annual  cycle,  the  diurnal  cycle,  and  the  first  harmonic  of 
the  diurnal  cycle,  were  obtained. 

A  complete  1-year  record  of  hourly  weather  data  was  then  syn¬ 
thesized,  using  the  mean  for  each  weather  variable  and  the  pair  of 
Fourier  coefficients  corresponding  to  the  annual,  the  diurnal,  and 
the  first  harmonic  of  the  diurnal  for  each  variable.  Thus,  at  the 
end  of  this  phase  of  the  analysis,  a  1-year  record  of  actual  weather 
data  for  1955  and  a  1-year  record  of  synthesized  weather  data  existed 
for  each  of  the  four  sites  selected. 

Building  Energy  Analysis  Test  Cases 

To  compare  the  effects  of  using  synthesized  and  actual  weather 
data,  a  hypothetical  building  was  designed  that  would  be  reasonably 
(but  not  overly)  sensitive  to  climatic  variations.  Also,  since  the 
performance  of  different  building  fan  systems  is  influenced  by  cli¬ 
mate,  even  when  satisfying  the  same  building  loads,  several  candidate 
building  fan  systems  were  designed  and  simulated.  Finally,  for  the 
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sake  of  completeness,  and  because  changes  in  hourly  demand  profiles 
can  cause  performance  to  vary,  a  boiler/chiller  plant  was  designed  to 
serve  the  building  fan  systems. 

The  hypothetical  building  consisted  of  24  exterior  office 

modules,  an  interior  computer  room,  and  three  interior  office 

modules.  The  exterior  modules  were  assumed  to  be  bounded  on  three 

sides  by  interior  building  partitions  and  were  assumed  to  have  an 

exterior  roof,  a  slab  floor  heavily  insulated  from  below,  and  an 

exterior  wall  of  various  constructions  and  with  varying  amounts  of 

glass.  Six  modules  were  assumed  to  be  on  each  of  the  north,  south, 

east,  and  west  sides  of  the  building.  Figure  5-2  illustrates  an 

office  module.  During  the  course  of  this  study,  three  different 

exterior  wall  constructions  were  considered:  one  was  all-glass, 

another  had  a  heavy  conductive  exterior  wall  with  about  16  percent 

glass,  and  the  third  had  a  heavily  insulated  exterior  wall  with  16 

percent  glass.  Table  5-1  gives  details  of  these  exterior  wall  con- 

2 

structions.  The  interior  office  module  was  assumed  to  be  72  m  in 
area,  square  in  shape,  and  in  contact  with  the  environment  only 
through  its  roof.  Three  identical  interior  office  modules  were  simu¬ 
lated.  The  remaining  module  represented  a  computer  room  and  was  dom¬ 
inated  by  interior  electric  load. 


Appropriate  occupancy  lighting  and  equipment  schedule  were 
selected  for  all  spaces,  and  a  control  strategy  designed  to  maintain 


Table  5-1 


Properties  of  Exterior  Walls,  Partitions,  Roofs,  and  Floors 


Layers 

Thickness 

(«) 

Dens  It v 
(kg/m3) 

Conduct lv ity 
(W/»  OK) 

Spec  111c 

Heat  (kJ/kg  ok) 

Solar 

Absorptivity 

Heavy 

Face  Brick 

.1020 

2002 

1.298 

.794 

.60 

Conductive 

Heavy  Concrete 

.1020 

7242 

1.731 

Hall 

Gypsum  Board 

.0123 

mu 

.640 

Heavy 

Wood 

.0254 

512 

1.379 

.115 

.  s 

1 nsulatcd 

Loose  Fill 

.0762 

5.4 

.038 

Wall 

Heavy  Concrete 

.1020 

224? 

1.731 

Glass  Hal  1 

Single-Pane 

.0063 

---{R  >  .OOB  m? 

<*/«>--- 

Heavy  Glass 

Partitions 

Concrete  8lock 

.2030 

689 

1.04 

.878 

Floor 

Insulat Ion 

. (R  • 

i«)  -  0». /U  I . 

Heavy  Concrete 

.  1020 

2242 

1.731 

.919 

Stone 

.0127 

881 

1.436 

1.674 

.55 

Membrane 

.0096 

1121 

.190 

1.674 

Insulat ion , 

.0509 

91 

.043 

.837 

Steel -Deck , 

.0015 

7688 

44.469 

.418 

Air-Space, 

-(R 

*  .  1 76  m2 

OK/W)--. 

Acoustic  Tile 

.0191 

480 

.061 

(he  space  t einpe  return  between  22.8  and  2‘j°C  was  assigned.  Healing 
and  cooling  were  assumed  to  be  available,  as  required,  throughout  the 
year. 


Five  separate  fan  systems  were  simulated,  with  each  designed  to 
serve  all  zones  in  the  building.  The  first  was  u  reheat  tan  system 
of  conventional  design.  Tin  second  was  a  reheat  fan  system  with  (1) 
cold-deck  reset  (i.e.,  the  cold  deck  was  reset  to  deliver  air  no 
colder  than  was  necessary  for  the  zone  requiring  the  most  cooling), 
and  (2)  an  enthalpy  economy  cycle  (any  time  the  enthalpy  of  the  out¬ 
door  air  was  less  than  the  enthalpy  of  the  return  air  for  the  system, 
outdoor  air  was  used  to  offset  all  or  part  of  the  cooling  require¬ 
ments).  The  third  system  was  a  conventional,  single-duct  variable 
air-volume  system  with  reheat.  The  fourth  system  was  a  variable 
air- volume  system  with  reheat,  cold-deck  reset,  and  an  enthalpy  econ¬ 
omy  cycle,  and  the  fifth  was  a  conventional,  four-pipe  fan-coil  sys¬ 
tem.  Figures  b-3  through  b-b  arc  schematics  of  these  systems.  Lach 
system  was  simulated  for  each  of  the  three  different  exterior  build¬ 
ing  treatments  previously  described.  Thus,  fifteen  1-year  fan  system 
simulations  were  performed  for  each  site  and  for  both  actual  and  syn¬ 
thesized  data.  Finally,  identical  modular  boiler/chiller  plants 
(three  chillers,  two  boilers)  were  simulated  for  each  simulated  fan 


system. 
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COIL  VAV  BOX 


volume  system 


HEATING  cooling 
COIL  COIL 


MOT  CtttllM 
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Analysis  of  Resul is 


We  now  assess  the  calculation  errors  in  expected  energy  use 
resulting  from  the  use  of  sinusoidal  weather  data  approximations. 

The  procedure  described  in  the  previous  section  provided  the  neces¬ 
sary  quantitative  data.  We  must  now  examine  this  data  to  determine 
whether  or  not  the  errors  are  acceptable,  lo  do  this,  we  will  review 
how  results  of  energy  analysis  calculations  are  used  (or  should  be 
used)  to  design  or  redesign  buildings  and  their  energy  systems  and 
determine  if  errors  resulting  from  the  use  of  synthesized  weather 
delta  affect  the  usefulness  of  the  results. 

Building  architecture  and  space  utilization  affect  building 
heating  and  cooling  leads.  Heating  and  cooling  loads  are  the  amount 
of  energy  which  must  be  added  to  or  removed  from  a  space  to  maintain 
it  within  the  specified  comfort  temperature  range.  A  design  objec¬ 
tive  is  to  minimize  these  loads  by  properly  selecting  the  building 
envelope  elements.  For  example,  we  might  expect  that  an  all-glass 
building  will  have  a  considerably  higher  cooling  load  on  a  hot  summer 
day  due  to  the  solar  transmission  through  the  glass  than  a  building 
with  a  modest  glass  area.  Similarly,  because  glass  is  not  a  good 
insulating  medium,  the  same  building  might  have  an  extremely  high 
h eating  requirement  on  a  cloudy  winter  day.  As  designers,  we  arc 
interested  in  both  the  heating  and  cooling  loads  within  a  building's 
zones  as  a  measure  of  the  building  envelope  performance.  If  we  are 
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to  use  a  s  i  ruiso  i  (Id  1  approx  iir.jl  ion  lot  climate,  simulation  results 
teem  tin  synthesized  dote  should  conform  in  results  obt  dined  by  Simu¬ 
la!  iny  the  zones  using  the  actual  weather  data;  in  particular,  the 
differences  between  different  building  architectural  and  construction 
options  should  be  the  same,  whether  we  use  synthesized  or  actual 
weather  data.  Figures  b-6  through  b-13  show  the  calculated  annual 
cooling  load  and  annual  heating  load  ior  the  exterior  zones  for  each 
of  the  four  climates  and  each  of  the  three  different  building  wall 
constructions  considered.  Results  of  the  use  of  actual  and  syn¬ 
thesized  weather  data  are  compared.  These  figures  provide  the  first 
evidence  that  we  are  on  relatively  firm  ground  when  we  treat  climate 
as  c  sum  of  sinusoidal  components.  Notice  in  particular  that,  except 
in  cases  where  the  differences  in  loads  are  very  small,  the  rank 
order  of  the  loads  for  each  option  remains  unchanged,  whether  we  use- 
actual  or  synttiesized  data.  For  example,  when  we  use  the  actual 
data,  the  annual  cooling  load  is  substantial ly  higher  for  an  all- 
glass  wall  option  than  tor  a  wel 1  - insu I ated  wall;  the  same  is  true 
when  we  use  the  synthesized  data.  Thus,  for  both  types  of  weather 
data,  the  difference  in  load  between  the  two  options  is  very  nearly 
the  same.  Lven  more  subtle  differences,  like  the  slight  increase  in 
calculated  annual  cooling  load  between  the  heavy-conductive  and 
heavy-insulat.ed  wall  options,  are  generally  preserved  when  syn¬ 
thesized  weather  data  is  used,  it  is  quite  clear  that  the  utility  of 
the  results  in  terms  of  assessing  the  performance  of  the  building 
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Figure  5-6.  Annual  cooling  load  for  exterior  zones,  Charleston,  SC 
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Figure  5-7.  Annual  heating  load  for  exterior  zones,  Charleston,  SC 
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Figure  5-8.  Annual  cooling  load  for  exterior  zones,  Madison,  WI 
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Figure  5-9.  Annual  heating  load  for  exterior  zones,  Madison,  WI 
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SYN  =  SYNTHESIZE 
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□  =  HEAVY  INSULATED  EXTERIOR 

Figure  5-10.  Annual  cooling  load  for  exterior  zones.  Fort  Worth,  TX. 
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Figure  5-11.  Annual  heating  load  for  exterior  zones,  Fort  Worth,  TX. 
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Figure  5-12.  Annual  cooling  load  for  exterior  zones,  Santa  Maria,  CA 
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Figure  5-13.  Annual  heating  load  for  exterior  zones,  Santa  Maria,  CA 
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envelope  is  the  same,  whether  we  use  synthesized  or  actual  data; 
there  is  only  a  small  error  in  absolute  value  if  we  use  the  syn¬ 
thesized  data. 

We  next  examine  the  fan  systems  meeting  heating  and  cooling 
loads  in  the  buildings  and  assess  whether  using  synthesized  weather 
data  will  cause  any  error  in  the  demand  that  the  fan  system  places  on 
the  central  plant  (i.e,  hot  and  chilled  water  demands).  Five 
separate  types  of  fan  systems  were  examined  for  each  architectural 
option.  Figures  5-14  through  5-21  summarize  the  results.  In  these 
figures,  system  1  is  the  conventional  reheat  system,  system  2  is  the 
reheat  system  with  cold-deck  reset  and  an  enthalpy  economy  cycle, 
system  3  is  the  conventional  variable  air  volume  system,  system  4  is 
the  variable  air-volume  system  with  cold-deck  reset  and  an  enthalpy 
economy  cycle,  and  system  5  is  the  four-pipe  fan-coil  system. 

We  must  now  assess,  for  a  given  building  envelope,  how  the  use 
\of  actual  vs.  synthesized  weather  data  will  affect  the  rank  order  of 
the  calculated  energy  consumption  for  each  fan  system  type.  An  exam¬ 
ination  of  Figures  5-14  through  5-21  clearly  indicates  that  there  is 
no  effect.  In  every  case  where  different  fan  systems  exhibit  sub¬ 
stantially  different  hot  and  chilled  water  demands,  the  variations 
are  roughly  the  same,  whether  we  use  synthesized  or  actual  data;  in 
fact,  we  see  better  agreement  in  the  absolute  values  for  annual  hot 
and  chilled  waste  demands  than  we  did  for  the  annual  heating  and 


ANNUAL  HEATING  COIL  DEMAND  ( GJ )  31  ANNUAL  COOLING  COIL  DEMAND  <GJ 


BLDG,  WITH  BLDG  WITH  HEAVY  BLOG.  WITH  HEAVY 

GLASS  EXTERIOR  CONDUCTIVE  EXTERIOR  INSULATED  EXTERIOR 


gure  5-14.  Annual  cooling  coil  demand,  Charleston,  SC. 
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Figure  5-15.  Annual  heating  coil  energy  demand.  Charleston,  SC 
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Figure  5-17.  Annual  heating  coil  demand,  Madison,  WI 
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Figure  5-18.  Annual  cooling  coil  demand.  Fort  Worth,  TX 


Figure  5-19.  Annual  heating  coil  ctemand,  Santa  Maria,  CA 
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5-20.  Annual  cooling  coil  demand,  Santa  Maria,  CA. 


GLASS  EXTERIOR  CONDUCTIVE  EXTERIOR  INSULATED  EXTERIOR 

Figure  5-21.  Annual  heating  coil  demand,  Santa  Maria,  CA 
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cooling  loads.  Thus,  using  either  synthesized  weather  data  or  actual 
weather  data,  building  designers  would  make  the  same  decision  in 
terms  of  the  design  of  the  building  fan  system. 

finally,  we  examine  the  performance  of  the  boiler/chiller  plant. 
In  the  context  of  this  simulation,  this  central  plant  is  the  final 
consumer  of  fuel  and  electric  power  and  is  the  point  where  we  can 
determine  the  cost/benefit  of  various  building  and  energy  system 
design  options.  Tables  5-2  through  5-5  present  overwhelming  evidence 
that  for  all  the  cases  studied,  use  of  either  actual  or  synthesized 
data  produces  the  same  rank  order  result.  After  the  synergistic 
effects  of  climate  on  the  building  loads,  fan  system  performance,  and 
central  plant  performance  have  all  had  an  opportunity  to  affect  the 
results,  the  absolute  value  for  energy  consumption  obtained  using 
either  actual  or  synthesized  data  are  nearly  identical.  Of  the  120 
comparisons  reported  under  the  "%  of  Actual"  heading  in  Tables  5-2 
through  5-5  (two  comparisons  per  plant  per  building  option  per  site), 
97  percent  of  the  synthesized  weather  data  results  are  within  5  per¬ 
cent  of  the  actual  weather  data  results.  Eighty-seven  percent  of  the 
results  agree  to  within  3-percent.  Notice  also  that  the  greatest 
percentage  disagreement  is  for  consumption  of  heating  fuel  at  sites 
like  Fort  Worth  and  Charleston,  where  very  little  fuel  is  required 
when  efficient  fan  systems  are  used.  In  these  instances,  a  rel  a- 
t.ively  small  error  in  absolute  value  is  a  relatively  large  fractional 
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Estimated  Total  Annual  Building  Energy  Consumption 
Charleston,  SC 


Buildiny  With  Glass  Exterior 


Total  Electric 

Fuel 

Plant 

No 

Actual 

SyillheSIZeU 

%  Of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

1311 

1311 

100.0 

1901 

19U3 

100.1 

2 

10/0 

1047 

97.6 

1284 

12/7 

99.5 

3 

1067 

106/ 

100.0 

700 

690 

98.6 

4 

9b6 

940 

98.3 

49b 

48/ 

98.4 

b 

939 

937 

99.8 

211 

203 

96.2 

building  With  Heavy 

Conductive  Extet ior 

lutal  tlictnc 

lucl 

Plant 

No 

Actual 

bynthcswid 

%  ot  Actual 

Ac  tua 1 

Synthesized 

t  ot  Actual 

1 

l  3U4 

1305 

100.0 

211/ 

2136 

100.9 

z 

1062 

1040 

97.9 

1471 

1478 

100.4 

3 

1041 

1042 

100.0 

910 

857 

94.1 

4 

919 

905 

98.5 

b86 

594 

101.3 

5 

89b 

094 

99.9 

196 

193 

98.5 

Bui  Idmg 

With  Heavy 

Insulated  Lxtenor 

Total  E let  trie 

Fuel 

Plant 

No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

l  of  Actual 

i 

1305 

130/ 

100.1 

1990 

2005 

100.7 

2 

1062 

1040 

97.9 

1342 

134? 

100.4 

3 

1042 

1042 

100.0 

721 

727 

100.8 

4 

918 

904 

98  5 

464 

467 

100.5 

5 

899 

899 

100.0 

78 

77 

96.3 

Not  cr 

— - 

— 

Plant 

1 

IS 

si  r  ving 

tilt  convent  tunal 

reheat  tan  system. 

Plant 

2 

is 

serving 

the  reheat  fan  system  vJtli  told-deck  reset  and  enthalpy  economy  eycle. 

Plant 

3 

IS 

si  r  v  ing 

the  conventional 

variable  air  volumi 

i  system. 

Plant 

4 

IS 

sir  v  i  ny 

the  variublt  air 

volume  system  with 

cold-deck 

reset  and  enthalpy  e-tonomy  < 

Plant 

5 

is 

si  r  v  i  ny 

the  conventional 

four  >plpc,  fan-coil  system. 

160 


Table  5-3 


Estimated  Total  Annual  Building  Energy  Consumption 
Madison,  WI 
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lot  al  l leclric 

Fuel 

Plant  No. 

Ac  t  uu 1 

by  nthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

1269 

1266 

99.8 

2560 
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2 
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940 

98.3 
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95.9 
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1044 
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98.9 
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98.8 

4 
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98.3 
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Actual 
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X  of  Actual 

1 

1266 

1264 

99.8 

2655 

2684 

101.1 

2 

932 

924 

99.1 

1888 

1913 

101.3 

3 

1016 

1014 

99.8 

1340 

1357 

101.2 

4 

1138 

829 

99.0 

1018 

1032 

101.3 

5 

885 

882 

99.6 

586 

586 

100.0 

Bu_i_l d i ng _ W i t li  Heavy 

Insulated  lxterior 

lotal  llectric 

F  uel 

Plant  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthes  r  zed 

X  of  Actual 

1 

1271 

1268 

99.8 

2340 

2350 

100.4 

2 

931 

924 

99.3 

1542 

1548 

100.4 

3 

1018 

1016 

99.8 

979 

984 

100.5 

4 

837 

830 

99.1 

686 

689 

100.4 

5 

889 

887 

99.7 

256 

253 

98.7 

1  -  * 
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Table  5-4 

Estimated  Total  Annual  Building  Energy  Consumption 
Fort  Worth,  TX 


Building  With  Glass  Exterior 


Plant  No. 

Actual 

Total  Electric 
Synthesized 

I  of  Actual 

Actual 

Fuel 

Synthesized 

i  of  Actual 

1 

1316 

1313 

99.8 

1856 

1826 

98.4 

2 

1086 

1063 

97.0 

1267 

1217 

96.1 

3 

1079 

1078 

99.9 

658 

632 

96.1 

4 

976 

966 

98.1 

478 

455 

95.1 

5 

969 

958 

99.8 

236 

204 

86.8 

Bui lding  With  Heavy 

Plant  No.  Actual 

Conductive  Exterior 

Total  Electric 

Synthesized  i  of  Actual 

Actual 

Fuel 

Synthesized 

t  of  Actual 

1 

1302 

1304 

100.1 

2037 

2060 

101.1 

2 

1069 

1041 

97.4 

1411 

1413 

100.1 

3 

1046 

1046 

100.0 

784 

789 

100.7 

4 

933 

914 

98.0 

542 

641 

99.8 

5 

909 

906 

99.7 

212 

194 

91.4 

Bui  lding  With  Heavy 

Insulated  Exterior 

Total  Electric 

Fuel 

Plant  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

l  of  Actual 

1 

1303 

1305 

100.1 

1935 

1953 

100.9 

2 

1068 

1041 

97.5 

1300 

1306 

100.5 

3 

1044 

1044 

100.0 

680 

683 

100.5 

4 

927 

909 

98.1 

429 

428 

99. 7 

5 

907 

908 

100.1 

87 

76 

87.3 

-  i 


Table  5-5 


Estimated  Total  Annual  Building  Energy  Consumption, 
Santa  Maria,  CA 


Bu i  I d i_ny  W i Id  01  ass  ExUrior 


Total  Electric 

fuel 

Plant  No. 

Actual 

Synt  lies  1  zed 

%  of  Actual 

Al  t  Ud  1 

Syntlu  sized 

X  of  Actual 

1 

lb4tl 

1942 

100.1 

3903 

3923 

100.6 

2 

899 

868 

96.9 

1629 

1564 

96.3 

3 

1049 

1049 

99.6 

913 

900 

98.6 

4 

820 

809 

98.6 

696 

673 

96.7 

5 

924 

921 

99.7 

285 

271 

96.0 

Buildinq  With  Heavy 

Conductive  Exterior 

Total  Electric 

f  ue  1 

Plant  No.  Actual 

Synthesized  %  of  Actual 

Actual 

Synthesized 

%  of  Ac 

1 

1287 

1287 

100.0 

2209 

2217 

100.3 

2 

861 

866 

99.3 

1424 

1427 

100.2 

3 

1033 

1032 

99.9 

841 

844 

100.3 

4 

792 

788 

99.5 

557 

559 

100.4 

6 

898 

897 

99.4 

102 

98 

96.4 

Bui Idi nj_Wi th  Heavy, 

,  Insulated  Ixltrior 

Total  Electric 

T  ue  1 

Plant  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

X  of  Actual 

1 

1287 

1287 

100.0 

2209 

2217 

100.3 

2 

861 

855 

99.3 

1424 

1427 

100.2 

3 

1033 

1032 

99.9 

841 

844 

100.3 

4 

792 

788 

99.5 

567 

559 

100.4 

5 

898 

897 

99.9 

102 

98 

96.4 

Use  of  synthesized  data  typically  produces  somewhat  smaller 
annual  heating  and  cooling  load  estimates  than  the  use  of  actual 
data;  however,  the  use  of  synthesized  data  never  affects  the  rank 
order  of  choices,  has  a  minimal  effect  on  the  difference  in  energy 
consumption  between  choices,  and,  in  the  final  analysis,  has  a  very 
small  effect  on  the  absolute  value  of  the  energy  use  calculated  using 
the  BLAST  energy  analysis  program.  In  fact,  the  use  of  synthesized 
data  need  not  produce  any  larger  differences  when  compared  to  actual 
data  than  the  use  of  two  different  years  of  actual  data.  Figures 
5-22  and  5-23  show  the  results  for  Madison,  W I ,  for  the  case  with 
heavy  insulated  walls  using  synthesized  data  from  1955  and  1960.  In 
many  cases,  the  difference  in  annual  load  estimates  made  with  the 
synthesized  data  compared  to  those  made  with  1955  data  are  less  than 
the  differences  between  the  1960  and  1955  weather  data. 

We  therefore  conclude  that  the  use  of  synthesized  climate  data 
composed  of  the  mean,  the  annual  sinusoidal  cycle,  the  diurnal 
sinusoidal  cycle,  and  the  first  harmonic  of  the  diurnal  cycle  for 
each  weather  parameter  provides  completely  adequate  weather  data  for- 
input  to  energy  analysis  tools.  Note  that  we  have  not  emphasized  the 
absolute  value  of  the  error  caused  by  using  synthesized  data  because 
it  is  only  in  very  specialized  applications  that  precise  absolute 
values  are  important  (i.e.,  establishing  conformance  with  building 
energy  performance  standards  --  a  legal  requirement).  Instead,  we 
have  focused  on  the  utility  of  the  results  for  comparing  system 


ANNUAL  HEATING  AND  COOLING  COIL  DEMANDS  (GJ) 


SYSTEM  1  SYSTEM  2  SYSTEM  3  SYSTEM  4  SYSTEM  5 


Figure  5-23.  Annual  heating  and  cooling  coil  demands  calculated  with 
synthesized  weather  data  and  with  weather  data  from 
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alternatives  --  the  principal  utility  of  energy  analysis  calcula¬ 
tions.  Surprisingly,  however,  there  are  only  small  errors  in  absolute 
value,  as  illustrated  by  Tables  5-2  through  5-5  and  by  the  extensive 
tabulations  in  Appendix  A,  which  show  the  percentage  errors  associ¬ 
ated  with  all  the  options  illustrated  in  Figures  5-6  through  5-21. 
These  errors  are  small,  even  though  we  have  not  included  higher-order 
harmonics  or  side  bands  in  the  weather  data  model. 

Having  concluded  that  the  use  of  a  sinusoidal  climate  approxima¬ 
tion  is  appropriate,  we  now  describe  the  load-calculating  scheme 
which  takes  full  advantage  of  such  an  approximation. 
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6  CALCULATING  HOURLY  LOADS  USING 
FREQUENCY  RESPONSE 

Introduction 

We  showed  in  Chapter  b  that  sinusoidal  weather  data  approxima¬ 
tion  produced  results  which  agreed  very  closely  with  those  of  actual 
data.  We  now  show  that  the  results  obtained  when  the  frequency 
response  approach  to  calculating  heat  conduction  through  multi¬ 
layered  slabs  is  used  agree  with  those  obtained  when  using  response 
factors  and  hourly  room  heat  balance  techniques. 

The  BLAST  program  was  used  to  calculate  the  hourly  loads  for 
several  typical  rooms  for  selected  days.  These  results  were  com¬ 
pared,  hour  by  hour,  with  results  calculated  using  a  sinusoidal  sol- 
air  temperature  as  the  outside  surface  temperature,  and  using  the 
frequency  response  conduction  model  described  in  Chapter  3.  Only 
heat  transfer  through  opaque  surfaces  was  considered,  since  in  the 
proposed  method,  solar  heat  gain,  and  internal  heat  gain  were  to  be 
accounted  for  using  existing  weighting  factor  techniques. 

Procedure 

The  BLAST  computer  program  calculates  room  heating  and  cooling 


loads  by  performing  a  detailed  heat  balance  on  both  the  inside  and 
outside  of  each  room  surface  and  on  the  room  air  volume.  For  a  room 


at  fixed-air  temperature,  the  heating  or  cooling  load  is  the  amount 
of  heat  which  must  be  added  or  removed  from  the  room  air  to  offset 
convective  heat  losses  or  gains.  While  convection  is  the  only 
mechanism  by  which  heat  can  be  transferred  to  the  room  air,  radiation 
and  conduction  must  also  be  considered  when  performing  the  heat  bal¬ 
ance  on  each  surface.  Radiation  and  conduction  affect  the  equili¬ 
brium  temperature  of  each  surface  and  thereby  affect  the  amount  of 
energy  transferred  by  convection  into  the  room  air. 

In  the  proposed  simplified  model,  heat  flux  through  opaque  walls 
and  roofs  is  calculated  using  an  approach  which  combines  the  effects 
of  radiation,  conduction,  and  convection.  For  each  exterior  wall  or 
roof,  we  assume  that  heat  transfer  takes  place  between  an  outside 
sol-air  temperature  and  a  fixed  inside  room  temperature.  The  sol-air 
temperature " i s  the  equilibrium  temperature  of  a  perfectly  insulated 
surface  of  specified  orientation  and  absorptivity  and  is  due  to  the 
combined  effects  of  incident  solar  radiation,  radiation  to  or  from 
the-  surroundings,  and  convection  to  or  from  the  ambient  air.  By 
using  the  sol-air  temperature  as  one  boundary  condition,  we  can  use  a 
combined  radiative-convective  outside  surface  heat  transfer  coeffi¬ 
cient  to  account  for  the  heat  transfer  by  both  mechanisms.  Simi¬ 
larly,  we  use  the  room  air  temperature  and  a  combined  inside  surface 
heat  transfer  coefficient  to  account  for  radiation  and  convection 
occurring  on  the  inside  surface.  We  first  calculate  the  thermal 
resistances  corresponding  to  these  inside  and  outside  heat  transfer 


coefficients;  these  are  used  with  the  thermal  resistance  and  capaci¬ 
tance  of  each  layer  in  the  wall  or  roof  to  calculate  the  frequency- 
dependent  transmittance  of  the  wall,  as  described  in  Chapter  3. 

Using  the  notation  of  Chapter  3,  the  real  part  of  the  frequency- 
dependent  transmittance  at  frequency  mis  UwS  and  the  imaginary  part 

is  Umc- 

We  now  show  by  simple  example  how  the  flux  is  calculated.  Sup¬ 
pose  the  sol-air  temperature  for  a  given  wall  or  roof  contains  a 
mean,  an  annual  sinusoidal  cycle,  a  diurnal  cycle,  and  the  first  har¬ 
monic  of  the  diurnal  cycle  (as  suggested  in  Chapter  4).  The  sol-air 
temperature  at  time  t  is  (use  Equation  4-1  for  comparison); 

TSUL_AiR(t)  =  aQ  +  2al  cos^t)  +  2bj  sin^t) 

+  2a365  cos(aJ365t}  +  2b365  sin^365t) 

+  2a730  cos(w730t)  +  2b73Q  sin(m730t)  [Eq  6-1] 

where : 


a 

a 

a 


(  is  the  mean  annual  sol-air  temperature 
and  b^  are  the  Fourier  coefficients  corresponding  to  annual 
sol-air  temperature  variation 

^  and  b^  are  the  Fourier  coefficients  corresponding  to  the 
diurnal  sol-air  temperature  variation 


&-J2 o  anc*  t*73Q  arc  the  Fourier  coefficients  corresponding  to  the 
first  harmonic  of  the  diurnal  sol-air  temperature  variation 
m.,  and  arc  the  frequencies  corresponding  to  the 

annual  cycle,  the  diurnal  cycle,  and  the  first  harmonic  of 
the  diurnal  cycle,  respectively. 

Since  even  the  heaviest  multilayered  walls  or  roofs  respond  in  steady 
state  to  the  annual  cycle,  we  can  lump  the  annual  mean  and  the  annual 
cycle  variation  into  a  local  mean  sol-air  temperature,  T  (t),  which 
i  s : 

Tm(t)  =  aQ  i  Zdl  cos(m1t)  +  2bj  sin^t)  Uq  6-2] 

The  flux  through  the  wall  is: 

q^  =  Uss^m^^  "  V  +  (U36bc  Zb365  +  U36bs  2a36b )cos(w365t ^ 

+  (U365s  2b36b  "  U36Sc  2a365^S1n^365t^ 

+  ^U73Cc  2b730  +  U7 30s  2d73C^c0s^73Ut^ 

+  ^U730s  2b730  '  U730c  2a730 ^ S1 n^W730t ^  b"3^ 

where : 

U36bs  and  U36bc  <irr‘  rca^  and  i,uayinary  parts,  respectively, 
of  the  transmittance  of  the  wall  at  the  diurnal  frequency 
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ll730s  dnc*  *^730c  are  t^ie  rea^  ancl  '’md9''nary  parts,  respectively, 

of  the  transmittance  of  the  wall  at  the  first  harmonic  of  the 
diurnal  frequency 

Uss  is  the  steady-state  transmittance  (U-value)  of  the  wall. 

Other  variables  are  as  defined  for  Equation  6-2. 

The  total  flux  through  the  wall  is: 


Q(t)  =  A  q(t) 


[Eq  6-4] 


where  A  is  the  area  of  the  wall. 


We  can  determine  the  total  flux,  QT(t),  due  to  heat  flow  through 
all  opaque  walls  and  roofs  by  applying  Equations  6-3  and  6-4  to  each 
wall  or  roof  and  summing  the  results.  Note  that  the  steady-state 
transmittance  and  the  components  of  transmittance  at  each  frequency 
are  different  for  each  different  wall  or  roof  construction  and  that 
the  Fourier  coefficients  for  sol-air  temperature  are  different  tor 
each  different  surface  orientation  and  absorptivity. 

To  determine  the  room  heating  or  cooling  load  due  to  heat  flow 
through  the  opaque  walls  and  roofs,  we  must  recognize  that  radiant 
heat  gain,  which  is  accounted  for  in  the  combined  inside  surface 
heat-transfer  coefficient,  is  really  a  delayed  room  load,  since  it 
must  be  absorbed  by  other  room  surfaces  before  it  can  be  convected 
into  the  room  air.  We  use  the  weighting  factor  approach  detailed  in 
the  ASHRAE  Handbook  of  Fundamentals  (ref  21)  to  account  for  this 


delay.  Using  this  procedure,  the  room  lead  at  time  t  caused  by  heat 
flow  through  opaque  walls  and  roofs  (or  by  conduction  but  not  by 
transmission  through  non-opaque  surfaces)  is: 
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L(0  =  VQ  Qj-(t)  +  Vj  Qt <  1-At )  -  WQ  L(t-At)  [Lq  6-5] 

where : 

L(t)  is  the  room  load  due  to  heat  flow  through  opaque  surfaces 
at  time  t 

L(t-At)  is  the  room  load  due  to  heat  flow  through  opaque  sur¬ 
faces  at  time  t-At  (the  previous  hour  if  At  is  1  hour) 

Qj(t)  is  the  total  heat  flux  through  opaque  surfaces  at  time  t 
Q-^(t-At)  is  the  total  flux  through  opaque  surfaces  at  time  t-At 
Vy,  Vj,  and  WQ  are  constants  which  depend  on  room  geometry  and 

construction.  Note  that  V,  =  1  +  W  -  V  . 

1  oo 

To  determine  if  any  substantial  error  results  from  using  Equa¬ 
tions  6-3  through  6-5  instead  of  the  detailed  heat  balance  embodied 
in  BLAST,  BLAST  was  executed  using  weather  data  which  yielded  sol-air 
temperatures  containing  only  a  mean,  a  diurnal  component,  and  the 
first  harmonic  of  the  diurnal.  The  room  loads  calculated  in  this  way 
were  compared  to  results  obtained  by  using  Equations  6-3  through  6-5. 
In  all  cases,  the  rooms  tor  which  loads  were  calculated  were  at  fixed 
inside  temperature  and  had  no  internal  or  solar  transmission  energy 
inputs.  The  same  room  geometry  and  wall  types  described  for  the 
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exterior  zones  in  Chapter  5  were  studied,  except  that  no  windows  were 
included  and  a  single  wood  sheet  was  used  in  lieu  of  the  all-glass 
exterior  as  representative  of  a  light,  conductive  exterior  wall. 

Resul ts 

Two  days  were  simulated:  one  was  a  cold  sunny  day  with  sol-air 
temperatures  on  a  south  vertical  exposure  (absorptivity  of  .6)  rang¬ 
ing  from  -14.5°C  to  15. 5°C ;  the  other  was  a  hot  sunny  day  with  south 
vertical  sol-air  temperatures  ranging  from  24°C  to  44.7°C.  Table  6-1 
provides  the  Fourier  coefficients  for  sol-air  temperature  for  these 
days.  These  sol -air  temperature  coefficients  were  used  along  with 
the  transmittances  shown  in  Table  6-2  to  calculate  the  heat  flow 
through  the  exterior  wall  and  roof  of  the  typical  exterior  zone 
described  in  Chapter  5.  Results  were  compared  to  BLAST  calculations 
made  with  appropriate  sinusoidally  varying  hourly  sol-air  temperature 
for  each  of  the  2  days  studied. 

Figures  6-1  through  6-7  compare  the  results.  For  cold  sunny 
days,  the  fluxes  and  loads  represent  outward  heat  flow  and  heating 
loads.  For  hot  sunny  days,  the  fluxes  are  heat  gains  and  the  loads 
are  cooling  loads.  In  each  figure,  the  curve  labeled  "Total  Flux"  is 
the  flux  calculated  using  Equations  6-3  and  6-4.  The  curve  labeled 
"Load,  Simplified  Model"  results  from  the  use  of  Equation  6-5.  The 


Table  6-1 


Fourier  Coefficients  for  Sol-Air  Temperature, 


South  Wal 1 , 

Absorptivity  - 

.6 

Mean 

a365 

b365 

a730 

b730 

Cold  Sunny 
Day 

1.3875 

-5.7417 

-4.4375 

1.0787 

.8237 

Hot  Sunny 

Day 

33.0417 

-3.9583 

-3.0592 

.7433 

.5650 

Table  6-2 

Frequency-Dependent  Transmittance  for  Exterior  Walls 

and  Roof 


Steady-State-., 
U-Value  (w/nrk) 

Diurnal 

1st  Harmonic 
Diurnal 

Section 

Real 

Imaginary 

Real 

Imaginary 

Heavy  Conductive  Wall 

2.9201 

.2065 

-1.5233 

-.4168 

-.6898 

Heavy  Insulated  Wal  1 

.4049 

.0110 

-.1551 

-.0339 

-.0715 

Wood  Sheet  Wal  1 

2.4646 

2.4420 

-.2684 

2.3760 

-.5264 

Root 

.4674 

.3455 

-.2594 

.1273 

-.3330 

POW 


0  2  4  6  8  10  12  16  18  20  22  24 

TIME  (HOURS) 


6-1.  Flux  and  loads  for  cold  sunny  day  and  a  heavy  conductive 
exterior  wall  using  precalculated  weighting  factors. 
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- TOTAL  flux 


. -  LOAD,  simplified  model 

—  —  blast  results 


TOTAL  OF  24  HOURLY  LOADS 


BLAST 

SIMPLIFIED 

MODEL 

28699 

28711 

■  »  -l  I  1  ±. -  1  I  1-  -  -1—  -1- 
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Flux  and  loads  for  cold  sunny  day  and  a  heavy 
conductive  exterior  wall. 


POWER  (WATTS) 


O  2  4  6  8  10  12  14  16  18  20  22  24 

TIME  (HOURS) 


Figure  6-3. 


Flux  and  loads  for  hot  sunny  day  and  a  heavy  conductive 
exterior  wal 1 . 
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Figure 


6-3.  Flux  and  loads  for  hot  sunny  day  and  a 
exterior  wal 1 . 


POWER  (WATTS) 


Figure  6-7.  Flux  and  loads  for  hot  sunny  day  and  a  wood 
exterior  walil . 
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"BLAST  Results"  are  the  heating  or  cooling  loads  calculated  by  the 
BLAST  energy  analysis  program. 

Figure  6-1  shows  Lite  results  for  a  cold  sunny  day  for  a  room 
with  a  heavy  conductive  wall.  While  the  total  of  the  hourly  loads 
calculated  using  the  simplified  model  is  nearly  the  same  as  the  BLAST 
total,  the  peak  load  is  higher  and  the  minimum  load  is  lower.  This 
discrepancy,  which  '  not  overly  severe,  was  determined  to  have 
resulted  from  the  use  of  pretabulated  weighting  factors  found  in  the 
ASHRAE  Handbook  of  Fundamental s .  However,  these  factors  do  not 
necessarily  apply  to  the  room  being  studied  (ref  22).  Adjustment  of 
these  weighting  factors  substantially  improved  the  results,  as  shown 
in  Figure  6-2.  Note  that  use  of  the  proposed  method  could  not  be 
justified  if  the  weighting  factors  were  adjusted  for  each  different 
day  type  analyzed;  however,  this  was  not  necessary,  as  shown  in  Fig¬ 
ure  6-3.  This  curve  is  for  the  same  zone  as  Figure  6-2  (heavy  con¬ 
ductive  exterior  wall),  but  is  for  a  hot  day.  The  same  weighting 
factors  were  used  to  calculate  loads  for  both  the  hot  and  cold  days. 

The  good  agreement  between  BLAST  and  the  simplified  method  cal¬ 
culations  shows  that  if  we  use  the  correct  weighting  factors,  the 
proposed  method  introduces  almost  no  error.  Figures  6-4  and  6-5  con¬ 
firm  this  result  for  a  heavy  insulated  exterior  wall;  Figures  6-6  and 
6-7  show  similar  excellent  agreement  when  the  exterior  wall  is  a  sin¬ 
gle  sheet  of  wood. 


Notice  that  we  have  tabulated  the  sum  of  the  hourly  loads  on 
each  figure  for  both  the  simplified  model  and  for  the  BLAST  calcula¬ 
tion.  These  results  are  remarkably  close.  In  the  worst  case,  they 
differ  by  only  2  percent. 

Conclusions 

The  preceding  analysis  shows  that  application  of  frequency 
response  techniques  produces  very  good  hourly  agreement  with  the 
response  factor  and  heat  balance  methods.  Recalling  from  Chapter  5 
that  the  annual  agreement  between  synthesized  and  actual  weather  data 
was  also  good,  we  may  conclude  that  the  application  of  frequency 
response  techniques  to  calculating  heating  and  cooling  loads  is  rea¬ 
sonably  precise  and  uncomplicated,  and  therefore  of  great  utility  for 
energy  analysis  programs. 
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7  CONCLUSIONS  AND  APPLICATIONS 

Conclusions 

We  now  summari zv  the  important  conclusions  which  can  be  drawn 
from  the  preceding  six  chapters. 

1.  In  Chapter  2,  we  showed  that  calculating  response  factors 
and  conduction  transfer  functions  is  tedious  and  computationally  dif¬ 
ficult.  Their  use  also  requires  a  substantial  amount  of  record  keep¬ 
ing  of  past  temperatures  and  fluxes  in  order  to  calculate  heat  con¬ 
duction  with  precision. 

2.  As  shown  in  Chapter  3,  calculating  the  frequency  response  of 
a  multi-layerr d  slab  (i.e.,  calculating  heat  flow  at  a  boundary  of 
fixed  temperature  caused  by  the  sinusoidal  variation  in  temperature 
on  the  opposite  boundary)  is  much  less  complicated  than  calculating 
response  factors  and  conduction  transfer  functions.  This  is  particu¬ 
larly  true  if  calculation  of  the  frequency  response  is  required  only 
at  3  few  selected  frequencies.  Attempts  to  simplify  the  calculation 
of  frequency  response  by  using  equivalent,  single-layer  slab  models 
do  not  appear  to  be  worth  t hr  error  and  risk  associated  with  such 
simpl if ications. 

3.  The  use  of  deterministic  plus  stochastic  modeling  tech¬ 
niques,  as  discussed  in  Chapter  4,  allowed  separation  of  the 
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deterministic  periodic  behavior  from  the  stochastic  behavior  in 
weather  data.  Important  periodicity  exists  in  the  weather  data  stu¬ 
died  only  at  the  annual  cycle,  the  diurnal  cycle,  and  harmonics  of 
the  diurnal  cycle.  These  are  the  only  frequencies  at  which  the  fre¬ 
quency  response  for  multi-layered  slabs  need  to  be  calculated. 

4.  The  use  of  deterministic  only  (sinusoidal)  models  for 
weather  data  accounts  for  a  large  fraction  of  the  variance  in  the 
data,  but  the  use  of  deterministic  plus  stochastic  models  accounts 
for  an  even  larger  fraction.  Even  though  we  have  chosen  to  ignore 
the  stochastic  autoregressive  behavior  in  the  proposed  simplified 
model  for  heating  and  cooling  load  calculations,  autoregressive 
behavior  is  definitely  a  characteristic  of  weather  data,  and  it  is 
appropriate  to  include  autoregressive  terms  in  weather  data  models 
used  for  other  purposes. 

5.  We  demonstrated  in  Chapter  5  that  the  use  of  synthesized 
weather  data  consisting  only  of  the  mean,  the  annual  cycle,  the  diur¬ 
nal  cycle,  and  the  first,  harmonic  of  the  diurnal  cycle  as  input  to 
the  BLAST  energy  analysis  program  produces  results  which  are  very 
nearly  the  same  as  results  produced  by  using  actual  weather  data  for 
a  wide  range  of  cases.  We  have  demonstrated  that  the  differences  in 
calculated  results  using  actual  weather  data  for  different  years  for 
a  particular  site  are  of  approximately  the  same  order  of  magnitude  as 
those  resulting  from  the  use  of  only  sinusoidal  weather  data 


j 


components  in  comparison  to  actual  weather  data.  We  may  conclude 
that  ignoring  the  stochastic  component  and  higher  harmonics  and  side 
bands  in  the  weather  data  produces  very  small  errors  m  calculated 
building  loads  and  energy  consumption. 

b.  the  use  of  frequency  response  techniques  to  calculate  con¬ 
duction  through  multi  layered  slabs,  combined  with  the  use  of  weight¬ 
ing  factors,  produces  reasonably  accurate  results,  hour  by  hour,  when 
compared  to  the  more  detailed  discrete  response  factor/heat  balance- 
methods  for  calculating  room  loads  (see  Chapter  6). 

7.  Using  the  results  of  Chapters  5  and  6,  we  may  conclude  that 
the  use  of  frequency  response  techniques  instead  of  the  more  detailed 
response  factor  and  heat  balance  methods  tor  calculating  air- 
conditioning  loads  produces  results  which  closely  agree  with  the  more 
detailed  methods.  The  procedures  are  so  much  simpler  that  they 
should  be  used  in  lieu  of  the  more  detailed  techniques. 

Applications 

The  utility  of  the  research  described  in  the  preceding  chapters 
is  that  it  establishes  the  technical  adequacy  of  a  load  calculation 
procedure  which  is  simple  enough  to  be  implemented  on  a  microcom¬ 
puter.  This  method  can  be  combined  with  existing  steady-state  tan 
system  and  boiler/chiller  plant  simulation  methods  (which  are  already 


relatively  simple)  to  produce  a  microprocessor-oriented  energy 
analysis  program. 

One  approach  to  such  a  program  is  shown  schematically  in  Figure 
7-1.  The  three  calculation  modules  shown  could  be  executed  sequen¬ 
tially  and  overwritten  in  the  microprocessor's  memory  with  the  next 
module  as  execution  progresses.  The  approach  shown  also  does  not 
require  storage  of  intermediate  results  which  will  minimize  storage 
requirements  and  I/O  operations. 

Two  sets  of  data  must  be  generated  separately  in  order  to  implt 
ment  a  microcomputer-based  energy  analysis  program.  First,  the  coef 
ficients  for  the  needed  weather  data  variables  must  be  generated 
using  the  methods  described  in  Chapter  4.  Second,  a  much  larger  set 
of  room  weighting  factors  than  are  shown  in  the  ASHRAF  Handbook  of 
Fundamental s  (ref  21)  should  be  tabulated.  Detailed  energy  analysis 
procedures  found  in  the  ELAST  program  should  be  used  to  calculate 
room  weighting  factors  for  a  wider  variety  of  room  geometries,  mess, 
and  glass  and  exterior  wall  areas. 

The  implementation  of  building  energy  analysis  procedures  on 
microcomputers  should  greatly  increase  the  number  of  engineers  who 
will  apply  these  procedures  to  improve  the  energy  efficiency  of  the 
buildings  and  energy  systems  they  design. 
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APPENDIX  A 

COMPARISON  OF  ENERGY  USE  PREDICTED  USING  ACTUAL 
AND  SYNTHESIZED  WEATHER  DATA 


Table  A-l 

Comparison  of  Energy  Use  Predicted  Using  Actual  and  Synthesized 
Weather  Data  for  Charleston,  SC,  and  for  a  Building  With 

Glass  Exterior 

Loaos 


Heating  Cooling 


Zone  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

19.26 

16.94 

88.0 

129.32 

126.06 

97.4 

2 

31.10 

28.88 

92.8 

126.42 

121.48 

96.1 

3 

32.60 

30.48 

93.5 

111.00 

106.14 

95.6 

4 

41.9/ 

40.68 

96.7 

84.18 

80.09 

95.1 

5 

0.0 

0.0 

100.0 

362.9 

362.9 

100.0 

t. 

3.78 

3.47 

91.9 

66.48 

55.34 

98.0 

fan  Systems 

Hot  Water 

Chil led 

Water 

System  No. 

A(  tua  1 

Synthesized 

l  of  Actual 

Actual 

Synthesized 

i  of  Actual 

1 

1327 

1324 

99.8 

2262 

2243 

99.2 

2 

866 

860 

99.3 

1185 

1061 

89.5 

3 

457 

463 

99.1 

1292 

1274 

98.6 

4 

315 

309 

98.1 

810 

735 

90.7 

5 

129 

120 

93.0 

881 

863 

98.0 

fioi  ler/Uil  1  ler  Plant 

Luc  1 

lotal 

1 1 ec 1 1  i  e 

I’lant  No. 

Actual 

Synthesized 

X  ot  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

1901 

1903 

100.1 

1311 

1311 

100.0 

2 

1284 

1277 

99.5 

1070 

1047 

97.8 

3 

700 

690 

98.6 

1067 

1067 

100 

4 

495 

487 

98.4 

956 

940 

98.3 

5 

211 

203 

96.2 

939 

937 

99.  H 

Table  A-2 


Comparison  of  Energy  Use  Predicted  Using  Actual  and  Synthesized 
Weather  Data  for  Charleston,  SC,  and  for  a  Building  With 
Heavy  Conductive  Exterior  Walls 

Looos 


Heat i nq 


Cool Ing 


Zone  No. 

Actua 1 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

X  of  Actual 

1 

22.40 

21.59 

96.4 

59.98 

66.69 

94.3 

2 

27.84 

26.89 

96.6 

61.52 

57.63 

93.7 

3 

28.60 

27.71 

96.9 

56.62 

53.54 

94.6 

4 

32.03 

31.30 

97.7 

50.08 

47.45 

94.7 

S 

0.0 

0.0 

100.0 

362.9 

362.9 

100.0 

6 

3.78 

3.47 

91.9 

66.48 

55.34 

98.0 

fan  Systems 

Hot  Water 

Chilled  Water 

System  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

1489 

1500 

100.8 

2223 

2215 

99.6 

2 

1017 

1025 

100.8 

1144 

1029 

89.9 

3 

562 

569 

101.3 

1166 

1157 

99.3 

4 

365 

370 

101.4 

655 

5  6/ 

89.6 

6 

115 

111 

96.5 

654 

639 

97.7 

Boiler/Chi  1 ler  Plants 

Fuel 

Plant  No.  Actual  Synthesized 

X  of  Actual 

Actual 

Total  tlectric 

Synthesized  I  of  Actual 

l 

2117 

2136 

100.9 

1304 

1305 

100.0 

2 

1471 

1478 

100.4 

1062 

1040 

97.9 

3 

910 

857 

94.1 

1041 

1042 

100.0 

4 

586 

594 

101.3 

919 

905 

98.5 

5 

196 

193 

98.5 

895 

894 

99.9 

Tabic  A-3 

Comparison  of  Energy  Use  Calculated  Using  Actual  and  Synthesized 
Weather  Data  for  Charleston,  SC,  and  for  a  Building 
With  Heavy,  Insulated  Exterior 

Loads 


Heating  Cool inq 


Zone  No. 

Actual 

Synthes i zed 

X  of  Actual 

Actual 

Synthesized 

X  of  Actual 

1 

6. SI 

5.95 

91.4 

67.82 

65.90 

97.2 

2 

9.35 

8.91 

95.3 

67.41 

65.32 

96.9 

3 

9.65 

9.33 

96.7 

63.46 

61.86 

97.5 

4 

11.33 

11.12 

98.16 

57.20 

55.92 

97.8 

5 

0.0 

0.0 

— 

362.88 

362.88 

100.0 

6 

3.78 

3.47 

91.9 

56.48 

55.33 

98.0 

Fan  Systems 

Hot  Water 

Chilled  Water 

System  No. 

Actual 

Synthesized 

X  of  Actual 

Actual 

Synthesized 

X  of  Actual 

1 

1396 

1403 

100.5 

2231 

2223 

99.6 

2 

926 

930 

100.4 

1143 

1028 

90.0 

3 

454 

437 

100.6 

1167 

1159 

99.3 

4 

272 

274 

100.8 

654 

587 

89.8 

5 

41 

39 

94.6 

683 

674 

98. 7 

Bo i  lei 

/Chi  1 ler  Plant 

Fuel 

Total 

L  lectric 

Plant 

No.  Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

X  of  Ai  t ual 

1 

1990 

2005 

100.7 

1305 

1307 

100.1 

2 

1342 

1347 

100.4 

1062 

1040 

97.9 

3 

721 

727 

100.8 

1042 

1042 

100.0 

4 

464 

467 

100.5 

918 

904 

98 . 5 

5 

78 

77 

98.3 

899 

899 

100.0 

4 
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Comparison  of  Fneryy  Use  Predicted  Usiruj  Actual  and  Synthesized 
Weather  Data  for  Madison,  Wl,  and  lor  a  Building  With 
(Hass  Exterior 


l  lid  CIS 


Hcdtjny  Cool_L0ii 


/out  No. 

At  l  uj  1 

Synthes i zed 

Of  Aitlldl 

At tu dl 

Synthes i zed 

% 

of  Actual 

i 

92.73 

81.13 

87.6 

f)M .  be 

86 .  /0 

87.9 

c 

117.81 

113.68 

96.6 

86.83 

76.36 

86.8 

i 

116.99 

112.79 

96.4 

86.13 

74.4 

87.4 

4 

13Z.00 

132.36 

100.3 

63.06 

4  7.99 

90.6 

b 

O.U 

0.0 

— 

363.4 

363.08 

99.9 

6 

16.98 

16.83 

99.1 

39.62 

38.60 

97.4 

fan  Systems 

Mot  Water 

Chil led  Water 

Systtm  No. 

Actual 

Synthes i 2ed 

X 

of  Actual 

Actual 

Synthes i zed 

X 

of  Actual 

1 

1838 

1846 

100.4 

1963 

1941 

99.4 

2 

139Z 

1329 

96.6 

706 

68/ 

83.2 

3 

9Z  / 

916 

98.7 

1 106 

1066 

96 . 4 

4 

740 

7Z0 

97.3 

463 

391 

84.4 

b 

4/7 

467 

96.8 

698 

6  74 

96 . 6 

Bo  i  1  cr/th i  i  1  er  Plant 

fuel 

loul 

tlectru 

Plant  No. 

Actual 

Synthesized 

X 

of  Actual 

Actual 

Syntht  s  i  zed 

X 

of  At  t  t>a  1 

1 

2660 

2678 

100.7 

1269 

1266 

2 

1991 

1909 

96.9 

966 

940 

♦if- . 

3 

1360 

1336 

98.8 

1044 

1032 

ii}.'  4  •. 

4 

1089 

1060 

97.3 

869 

866 

Oy  .  • 

s 

704 

676 

96.1 

917 

91  1 

uu .  4 

Tabl  e 
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Comparison 

of  Energy  Use  Predicted 

Using  Actual  and  Synthesized 

Weather 

Data  for  Madison,  WI , 

and 

for  a 

Building  With  Heavy 

Conductive 

Exterior 

Loads 

Heating 

Cool i ng 

Zone  No.  Actual 

Synthesized 

l  of  Actual 

Actual 

Synthesized 

%  of  Actual 

I  84.32 

82.68 

96.1 

38.82 

33.05 

85.1 

2  93.61 

93.53 

100.0 

36.96 

32.90 

96.9 

3  93.35 

93.16 

99.8 

36.71 

32.31 

88.0 

4  99.06 

99.66 

100.6 

29.03 

26.46 

91.1 

5  0.0 

0.0 

— 

353.41 

363.09 

99.9 

6  16.98 

16.82 

99.0 

39.62 

38.6 

97.4 

Fan  Systems 

Hot  Water 

Chilled  Water 

System  No.  Actual 

Synthesi zed 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

936 

1928 

99.6 

536 

498 

492.9 

992 

980 

98.8 

313 

276 

88.2 

537 

620 

96.9 

Boiltr/Chiller  Plant 

Fuel 

Total 

tlectric 

Plant  No.  Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

Table  A-6 


Comparison  of  Cnergy  Use  Predicted  Using  Actual  and  Synthesized 
Weather  Data  for  Kadi  son,  WI,  and  lot  a  Building  With  Heavy, 

Insulated  l.xterior 

Loads 


Med  1 1  riy  (  ooluig 


Zoiu*  No. 

AtllJol 

Synt  hesi  z«  d 

X 

of  At  t  u<j  J 

A*  Mi.il 

Synt  lies  i  zed 

X 

of  Am. u ril 

i 

Z8.es 

Z7.37 

94  .9 

4  i .  60 

46.6 

9b.r 

? 

33.62 

33 . 39 

98.7 

46.04 

42.72 

44.9 

3 

33. 70 

32.98 

97.9 

44.66 

42.23 

94 .  < 

4 

3b. 33 

36. 36 

100.0 

17.41 

31 .94 

% .  fj 

S 

0.0 

0.0 

— 

363.4 

363.08 

99 . 9 

b 

16.97 

16.82 

99 . 1 

:<9.w 

94  / .  4 

fan  Systems, 

Hut  Water 

tin  1  ltd  Water 

System  No. 

Actual 

Synthesized 

X 

ol  Actual 

A?  tudl 

Synthes i zed 

Cf> 

o  1  Ac  t 1 ,  ci  1 

1 

164  7 

1667 

100.6 

I960 

1449 

"9.4 

2 

10/8 

1084 

100.6 

1,  U, 

49K 

"7.1 

J 

663 

666 

100.4 

1002 

990 

W.* 

4 

43b 

4  36 

100.1 

HZ 

2/4 

1 .4 . .( 

t, 

ISO 

146 

96.4 

672 

66  1 

OH .  4 

Bo l  It  r/Lhi  1  ler  Plant 

Fuel 

To!  a  1 

L  let  It  ic 

Plant  No. 

Actual 

Synthesized 

% 

of  Actual 

Ac  tua 1 

Synthesized 

% 

Of  fit.  1  1.10  1 

Table  A-7 


Comparison  of  Energy  Use  Predicted  Using  Actual  and  Synthesized 
Weather  Data  for  Fort  Worth,  TX,  and  for  a  Building 
With  Glass  Exterior 


Loads 


Heat i ng 


Cool  mg 


Zont  No. 

Actual 

Synthes  1  zed 

%  of  Actual 

Actual 

Synthesized 

T  of  Actual 

1 

22.96 

15.70 

68.4 

153.04 

151.55 

99.0 

2 

35.50 

29.55 

83.2 

154.29 

146.6 

95.0 

3 

36.46 

30.65 

84.1 

139.95 

132.8 

94.7 

4 

46.89 

42.66 

91.0 

99.50 

93.85 

94.3 

5 

0.0 

0.0 

— 

364.68 

364.32 

99.9 

6 

4.38 

3.53 

80.6 

61.68 

60.14 

97.2 

Fan  Systems 

Hot  Water 

Chilled  Water 

System  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

t  of  Actual 

1 

1296 

1270 

98.0 

2299 

2255 

98.1 

2 

849 

814 

95.8 

1266 

1091 

86.0 

3 

435 

417 

95.9 

1351 

1322 

97.9 

4 

312 

292 

93.5 

909 

803 

88.3 

5 

146 

122 

83.6 

987 

962 

97.5 

Boi ler/Chi 1 1 er  Plant 

Fuel 

Total 

Electric 

Plant  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

1856 

1826 

98.4 

1316 

1313 

99.8 

2 

1267 

1217 

96.1 

1086 

1053 

97.0 

3 

658 

632 

96.1 

1079 

1078 

99.9 

4 

4/H 

455 

95.1 

9/*, 

956 

98. 1 

ij 

2.15 

204 

86. 7 

•IVI 

958 

99,8 

Table  A-8 


Comparison  of  Energy  Use  Predicted  Using  Actual  and  Synthesized 
Weather  Data  for  Fort  Worth,  IX,  and  for  a  Building 
With  Heavy  Conductive  Exterior 


Loads 


Heat  i  ng  "  Coo 


/unt  No. 

A(  tu.i  1 

Synthesized 

%  of  Actual 

Act  ual 

Synthesized 

5  of  Actual 

1 

25.12 

20.  OH 

83.1 

/6.8b 

.'1.28 

92 .  / 

Z 

30.6/ 

2/. 00 

88.0 

80.11 

/2.90 

91 .0 

3 

31.30 

2/. 6b 

88.3 

/4.5H 

68.69 

92 . 1 

4 

34.9/ 

28.12 

80.3 

63.33 

by.// 

92.8 

5 

0.0 

0.0 

— 

364./ 

364.3 

99.9 

6 

4.38 

3.53 

80.6 

61 .88 

60.14 

97.2 

Can  Systems 

Hot  Water 

Chi  1  led  Water 

System  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthes i zed 

5t  Of  Act  ual 

1 

14  34 

1446 

100.9 

2218 

2200 

99.18 

2 

96b 

9/6 

100.8 

11/8 

1031 

87.5 

3 

blh 

51/ 

100.2 

1192 

1170 

98.1 

4 

345 

340 

98.  / 

721 

624 

86.5 

S 

1  26 

110 

>■ ) .  3 

730 

704 

96.4 

lio  l  U  r  / O  i  l  Ur  Plant 

1  UL  1 

Total 

tlectric 

21  ant  No. 

Art  tin  1 

Synt  host  Zed 

J,  of  Ac  tual 

Act  ual 

Synthesized 

t  ot  Actual 

1 

203/ 

2060 

101.1 

1302 

1304 

100.1 

2 

1411 

1413 

100.1 

1069 

1041 

Cl  7  4 

3 

/84 

789 

100.7 

1046 

1046 

100.0 

4 

542 

541 

99.8 

933 

914 

98.0 

b 

212 

194 

91.4 

909 

906 

99.7 

Table  A-9 


Comparison  oi  Energy  Use  Calculated  Using  Actual  and  Synthesized 
Weather  Data  for  Fort  Worth,  TX,  and  for  a  Building  With  Heavy, 

Insulated  Exterior 

Loads 


Heating  Cool ing 


Zone  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthes l zed 

l  oi  Actual 

1 

Z .  7  0 

5.61 

72.9 

77.97 

75.32 

96.6 

2 

10.57 

8.99 

85.1 

78.82 

74.82 

94.9 

3 

10.70 

9.29 

.  86.8 

74.24 

71.26 

96.0 

4 

12.53 

11.43 

91.2 

64.11 

62.01 

96.7 

S 

0.0 

0.0 

— 

364.68 

364. 32 

99.9 

6 

4.38 

3.53 

80.6 

61.88 

60.14 

97.2 

Fan  Systems 

Hot  Water 

Chi  1  led  Water 

System  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

X  of  Actual 

1 

1360 

1370 

100.7 

2224 

2207 

99.2 

2 

895 

900 

100.5 

1175 

1028 

87.5 

3 

428 

429 

100.2 

1181 

1162 

98.4 

4 

254 

252 

99.3 

697 

604 

86.6 

5 

46 

39 

84.4 

730 

716 

98. 0 

Boylei 

/Chiller  Pl_ant 

1  1IC  1 

total 

Llectnc 

Plant 

No.  Actual 

Synthesized 

l  oi  Actual 

Actual 

Synthesized 

X  oi  Actual 

1 

1935 

1953 

100.9 

1303 

1305 

100.1 

2 

1300 

1306 

100.5 

1068 

1041 

97.5 

3 

680 

683 

100.5 

1044 

1044 

100.0 

4 

429 

428 

99.7 

927 

909 

98.1 

5 

87 

76 

87.3 

907 

908 

100.1 

Tabic  A- 1C 


Comparison  of  Energy  Use  Predicted  Using  Actual  and  Synthesized 
Weather  Data  for  Santa  Maria,  CA,  and  for  a  Building  With 

Class  Exterior 

ludds 


Heajting  too]  iny 


/onr  No. 

Act  uu  1 

Synthesized 

%  oi  Actual 

A.  tu.il 

Syn!  in  1  Z<  .1 

%  of  At  Ti.<*  1 

i 

2K.oi 

23.39 

83.6 

106.1/ 

102.97 

0  /  p  u 

2 

41.97 

39.  HI 

94.8 

70.70 

71.82 

U.'j  ,  f 

i 

32.SU 

28.  SI 

87.7 

92.92 

82.04 

Kl.j 

4 

56.27 

56.7 

100.7 

36.18 

33. IS 

91  .(■ 

s 

0.0 

O.U 

— 

3S8.12 

JS7.91 

<HJ  .  U 

6 

S.21 

4.98 

95 .6 

42.74 

41.87 

91.! 

Fiiit  Systems 

Hot  Water 

Chi  1  led  Water 

System  No. 

Actual 

Synthes i zed 

%  of  Actual 

Actual 

Synthesi zed 

1  of  Actual 

1 

2S4 1 

2S49 

100.3 

2941 

293S 

on .  j: 

2 

1136 

1090 

96.0 

364 

2S2 

bu.? 

3 

610 

S99 

98.3 

1082 

1062 

9K2 

4 

443 

426 

96.2 

214 

1 62 

;n.h 

S 

164 

l  S3 

93. S 

702 

688 

w.l 

Boi 1 er/lhi tier  Plant 

Fuel 

lotal 

Electric 

Plant  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

3  SO  3 

3S23 

1  (JO .  6 

IMO 

1 S42 

KK.l 

/ 

1  025 

1 564 

% .  3 

p.nu 

068 

ut  . 1 

j 

91  ( 

900 

W .  (, 

KM'i 

1045 

1 

1.96 

6/  ! 

%.  / 

0.09 

i 

s 

285 

2/1 

921 

Tabic  A- 11 
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Comparison  of  Energy  Use  Predicted  Using  Actual  and  Synthesized 
Weather  Data  for  Santa  Maria,  CA,  and  for  a  Building  With  Heavy, 
Conductive  Exterior  Walls 

Luads 


Heat  ing  Cool  i  ng 


/om  No. 

Actud 1 

Synthes i zod 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

i 

33.su 

31.69 

94.6 

29-/4 

27.29 

91.7 

2 

40.89 

40.48 

99.0 

25.70 

23.6/ 

92.1 

3 

36.63 

35.81 

97.8 

27.28 

23.82 

87.3 

4 

48.25 

48.49 

100.5 

17.54 

16.12 

91.9 

5 

0.0 

O.U 

— 

358.12 

357.91 

99.9 

6 

5.21 

4.98 

95. b 

42.74 

41.87 

97.7 

Fan  Systems 

Hot  Water 

Chilled  Water 

System  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

1720 

1727 

100.4 

1972 

1969 

99.8 

2 

1173 

1179 

100.5 

183 

138 

75.6 

3 

725 

730 

100.7 

940 

991 

100.1 

4 

603 

505 

100.5 

98 

69.7 

71-1 

5 

165 

161 

97.8 

604 

493 

97.9 

Boi ler/Chi 1  ter  Plant 

fuel 

Total 

Electric 

Plant  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

%  of  Actual 

1 

2449 

2460 

100.4 

1284 

1284 

100.0 

2 

1654 

1656 

100.1 

M(>1 

H6S 

99.3 

J 

1081 

1 089 

1(10. 11 

Itl.tl 

10.11 

100.0 

4 

/90 

7% 

100.  / 

/92 

/M'l 

99 .0 

5 

297 

294 

98.9 

8.84 

1182 

99.8 

Table  A- 12! 


Comparison  of  Energy  Use  PredicLed  Using  Actual  and  Synthesized 
Weather  Data  for  Santa  Maria,  CA,  195b,  and  for  a  Building 
With  Heavy,  Insulated,  Exterior  Walls 


Loads 


Heating 


Loot  i ng 


l ont  No. 

Actual 

Synt  Iiesi2e(l 

t  of  Actual 

Actual 

Synt  hi  s  i  zed 

%  of  /W.tual 

1 

8.40 

7.32 

87.1 

50.80 

49.65 

9/. 7 

2 

11.96 

11.64 

97.3 

45.53 

44.32 

97.3 

3 

10.49 

9.72 

92.6 

49.14 

46.65 

94.9 

4 

14.89 

15.00 

100.7 

36.94 

36.03 

94.9 

5 

0.0 

0.0 

— 

358.12 

6 

5.21 

4.98 

95.5 

42.74 

41  .87 

97.7 

fan  Systems 

Hot  Water 

(.hi  1  led  Water 

System  No. 

Actual 

Synthesized 

%  of  Actual 

Actual 

Synthesized 

X  of  Actual 

1 

1541 

15A6 

100.3 

1987 

1984 

99.8 

2 

997 

1000 

100.3 

183 

138 

75.6 

3 

540 

542 

100.5 

1005 

1001 

99.6 

4 

335 

336 

100.  J 

105 

76 

72.6 

5 

51 

49 

95.4 

590 

581 

98.5 

Boi lei/thi 1 ler  Plant 

Fuel 

Total 

Llectric 

Plant  No. 

Actual 

Synthes i zed 

*  of  Actual 

At  tual 

Synthesized 

X  of  Actual 

1 

2209 

2217 

100.3 

1287 

1287 

100.0 

2 

1414 

1427 

100.2 

861 

855 

99.3 

3 

841 

844 

100.3 

1033 

1032 

99.9 

4 

557 

559 

100.4 

792 

788 

99.5 

5 

102 

98 

96.4 

898 

897 

99.9 
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VITA 

Douglas  C.  Hittle  was  born  on  April  IP,  194/, in  Fort  Collins, 
Colorado.  Fie  earned  a  Bachelor  of  Science  degree  in  Mechanical 
Engineering  in  1969  and  a  Master  of  Science  degree  in  Environmental 
Engineering  in  1975,  both  at  the  University  of  Illinois,  Urbana- 
Champaign. 

Mr.  Hittle  has  engaged  in  a  full-time  engineering  career  since 
June  1969,  when  he  joined  the  Base  Civil  Engineer's  staff  at  Chanute 
Air  Force  Base,  Rantoul  ,  Illinois,  as  a  mechanical  engineer  having 
design,  inspection,  and  trouble-shooting  responsibilities  for  large 
air-conditioning  and  other  mechanical  systems.  In  June  1973,  he 
joined  the  U.S.  Army  Construction  Engineering  Research  Laboratory 
(CERL),  Champaign,  Illinois,  where  he  now  serves  as  a  team  leader  in 
the  Energy  Systems  Division.  At  CERL,  Mr.  Hittle  has  been  responsi¬ 
ble  for  providing  technical  guidance  to  research  teams  performing 
energy-related  research  and  development  including  the  development  of 
the  Building  Loads  Analysis  and  System  Thermodynamics  (BLAST)  com¬ 
puter  program  for  simulating  energy  use  in  buildings  and  the  develop¬ 
ment  of  simplified  procedures  for  predicting  the  performance  of  solar 
energy  systems. 

Mr.  Hittle  was  named  the  CERL  Researcher  of  the  Year  in  1976, 


received  the  Department  of  the  Army  Research  and  Development 


Achievement  Award  in  1977,  and  was  awarded  a  Year  of  Advanced  Study 
in  1979.  He  is  a  member  of  Sigma  Xi  research  honorary  society. 

Mr.  Hittle  is  also  active  in  the  American  Society  of  Heating, 
Refrigeration,  and  Air-Conditioning  Engineers,  where  he  serves  on  th 
Technical  Committee  on  Energy  Calculation,  the  Subcommittee  on  Sim¬ 
plified  Energy  Cal cul at  ions.  Panel  10  of  the  Standards  Committee  for 
Energy  Conservation  in  New  Buildings,  and  as  chairman  of  the  Subcom¬ 
mittee  on  System  Simulation.  He  is  also  a  member  ol  the  subgroup  on 
Heating  and  Cl imat ization  of  the  International  Council  for  Building 
Research  Studies  and  Documentation  (C1B). 

Mr.  Hittle's  publications  include: 
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to  Buildings,  May  1978. 
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Locations,  ASRAE  Trans.,  Vol.  83,  Part  I,  1977. 
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tion  in  Residential  and  Commercial  Buildings,  ASHRAE  Trans.,  Vol .  83, 
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% 

Technical  Manuscript  E-169,  " Calculating  Building  Heating  and  Cooling  Loads 
Using  the  Frequency  Response-of  Multilayered  Slabs",  by  D.  C.  Hittle 
(February,  1981)  ADA097597. 


Page  24;  [Eq  2-16] 

Replace:  B  =  (sinh  (i  /s/a  ))  T2(s) 


coshU  /s/a) 


k 


wi  th : 


B  »  (sinh  (i  /~s7a))  T2(s) 


coshU  /s/a) 


k 


92(s) 


Page  32:  fi rst  line: 

Replace:  R„C 
n  pn 


with: 


R_C 
n  n 


Pace  34:  figure  2-3 
Replace:  labels  of  dashed  lines: 


r'iur«  2-j.  "•»»  -r’tnqulir  j«i  tt  n  sum  of  :nt  :hrtt  r mios. 
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Page  37:  [Eq  2-42] 

Replace: 

-  /sOT  est 
r  a  —  ,  1  1 

0  r/s^i C-j 


with: 


-  /sR-j  C1  eSt 
n  ^Rls2  sinh 


r  = 


s*- 


-n2ir2 
R1 C! 


Note  lower  case  s 


/-n2ir2 

g-tn^/RjCj 

TnV 

cosh  /-n2m2 

1  ¥7 

2/-n2n2 

•/-n2ir2 

e-tn2x2/R1C1 

.  nV 

cosh  /-n2m2 

1  ¥i 

2/-n2n2 

Page  37:  just  after  [Eq  2-42]  the  phrase: 
Replace  the  phrase: 

/ 

Now,  cosh  /-A*  *  cos  (jnir)  *  (-l)n 


Wi’n"  Now,  cosh  f-n2*2  »  cos  (rnr)  *  (-1) 

Page  39:  just  after  [Eq  2-49] 

Replace  the  phrase: 

Similarly,  for  any  incremental  time,  m6,  where  m  is  greater  than  or 
equal  to  3: 

with:  Similarly,  for  any  incremental  time,  m  5  ,  where  m  is  greater  than 

equal  to  3:* 


and  add  footnote: 

Mt  is  actually  valid  for  m  *  2  also,  but  the  convergence  is  too  slow  to 
make  it  useful;  the  remainder,  which  is  dropped,  is  proportional  to  1  . 
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Page  47: 
Replace: 

with: 

Page  49: 
Replace: 


with: 


Page  51 : 
Replace: 

wi  th : 

Page  51_ 
Replace 


[Eg  2-66] 

q,(2S)  *  I,X2  -  T0X2 


qi(24)  -  t,xz  -  y2 

footnote  at  bottom 


P(s)  * 

6S 

P(s)  *  (1-2e“s6)/6s2 

P(s)  -  (l-2e‘s<S  ♦  e‘2s5)/Ss2 


0  <_  t  <  5 
6  <  t  <_  26 

t  >  26 


P(s)  =  (l-2e”s<5  +  e”2s<5)/ 5s2 


first  line 

This  is  the  positive  X-directed. .. . 
This  is  the  positive  x-directed.... 


[Eg  2-76] 


8  ( s ) 


*1 

+  —  1  ■■  cosh  /s"R9C» 


'T^T1CHP* 


mu.  wpiMpppi 
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+ -  cosh  /sR-C-  sinh  /s  R. C, 

ZsR^  22  11 

(Insert  s) 

Page  53:  just  above  [Eq  2-80] 

Replace:  The  y  response  . 

with:  The  Y  response  . 


with: 

B(s)  = 


Page  55: 
Replace: 


top  line 


_d_ 

ds 


'dA(s) 

ds 

dc(s) 
L  ds 


Mil 

ds 

Mil 

ds 


wi  th : 


_d_ 

ds 


■■ 

r«iii 

Mil' 

p 

_ 

ds 

ds 

s 

_ 

Mil 

L  ds 

Mil 
ds  J 

Page  56: 
Replace: 


[Eq  2-87] 
qo,t  =  T  i  y 


m 


Note  upper  case  C 


wi  th: 

Page  57: 
Replace: 

wi  th : 


Page  57: 
Replace: 

wi  th : 

Page  58: 
Replace: 

with.: 

Page  59: 
Replace: 


IP,  pJMW  -  - 
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i 


’o,t  =  T.  I  y 


m»l 


m 


mid-page  paragraph:  (3rd  line  of) 


Recall  that  n  . 

on  the  . 

increases  . Rj  =  .4m2.  °K/w  and 

Recall  that  n  . 

on  the  . . . 

increases  . R.  =  .4  „  0K/W  and 


last  line 


6 

-126.17 

x  10’5 

-4.5423 

.0149 

-9500.7 

6 

-126.17 

x  10"5 

-4.5423 

.0106 

-9500.7 

5th  line 

contribution  caused  .  =  ^1  1  +  92A2* 

contribution  caused  .. .  ^  a  ^  +  ^2^2* 


4th  line  from  the  bottom 


Notice  that 


Tq,  t-m+1 
T0,  t-m'+l 


with: 


Notice  that 


Page 

i  61: 

5th  line 

Replace: 

factors  until  . 

with: 

factors  unti 1  . 

Pag< 

»  64: 

[Eq  2-112] 

Replace: 

m" 

<,  t  *  l  . 

m*l 

.  *  (VVqi,t-rVi,t-2 

+  (X1+X2)q.  #t.1-X1X2  qi>t_2 

•  'Xlqo,t-2 

•  "X?  X2qo,t-2 


with: 


with: 
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Page  68: 
Replace: 

with: 

Page  69: 
Replace: 


last  line 

the  test  of  Equation  2-125  is  applied  for  successively  increasing 
the  test  of  Equation  2-126  is  applied  for  successively  increasing 

second  line 

made  for  k  <  M.  As  . Equation  2-125 . . . 


with: 


made  for  k  <  M.  As 


Equation  2-126 


